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Part I: Statement of the Problem and solution technique 


1. INTRODUCTION 


1.1. Background 

Hie distinct feature of a Stirling engine as compared u roost other power productng 
devices is its ability to be driven by virtually any heat source such as solar energy or 
combustable refuse, among others. This leads to a number of promising applications: Hie 
Stirling engine as a prospective power source for future space missions or as the rice husk 
driven motor for agricultural machinety. In addition, a Stirling engine is silent, has low 
emission levels if powered by a combustion process and is energy efficient. 

A crucial point for further development of Stirling engines is the optimisation of its heat 
exchangers whtch operate under oscillatory flow conditions. Heat exchanger optimiiation 
depends on the ability to accurately simulate the fluid flow and heat transfer behavior. It has 
been found that the most important thermodynamic uncertainties in tire Stirling engine 
designs for space power are in the heat transfer between gas and metal in all engine 
components and in the pressure drop across the heal exchanger components. So far, 
performance codes cannot predict the power output of a Stirling engine reasonably enough 
if used for a wide variety of engines Thus, there is a strong need for belter performance 
codes. How ever, a performance code is usually not concerned with the details of the flow . 
This information must be provided externally. While for laminar oscillating flow analytical 
relationships exist, there has been hardly any infomtation about transitional and turbulent 
oscillating flow w hich could be introduced into the performance codes. 
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In 1986 a survey by Seume and Simon (1986a) revealed that most Stirling engine heat 
exchangers operate in the transitional and turbulent regime. Consequently, research has 
since focussed on the unresolved issue of transitional and turbulent oscillating flow and 
heat transfer. Since 1988, the University erf Minnesota oscillating flow test facility has 
obtained experimental data about transitional and turbulent oscillating flow. However, since 
the experiments in this field are extremely difficult, lengthy and expensive, it will be 
advantageous to numerically simulate the flow and heat transfer accurately from first 
principles. With a simulation program, one can enhance the understanding of the oscillating 
flow phenomena in general. Also, a simulation program is useful in guiding the experiment 
in some areas. Once tested for its validity, many more operating points can be “probed" 
with such a simulation program than with any experimental set-up. Finally, this program 
can generate the input data needed for a performance code, as mentioned above. 

It is the purpose of this research to contribute to the development of a realistic 
simulation program, thereby adding to the basic knowledge and understanding of turbulent 
oscillating flow and heat transfer and to the further the development of Stirling engines. 

1.2. Literature Survey of Oscillating Flow Research 

1.2.1. Scope and limitations of this Survey 

The objective of this survey is to answer the following questions: 

On oscillatory pipe flow, ... 

a) ... what experimental data is available? 

b) ... what numerical results are reported and how do they compare with 

experiments and theory? 

c) ... what are the open questions for numerical analysis in this field? 
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Currently there is no consensus in the literature about the nomenclature of periodic 
unsteady flows. Expressions like pulsating, pulsatile, oscillating or oscillatory are used 
synonymously. In this work, we are primarily concerned with information about 
periodically reversing flow with zero mean velocity. In the following, we shall denote this 
situation by "oscillatory" or "oscillating" flow only. Tbe expressions "pulsatile" or 
"pulsating" flow will be used for flow situations with a net mean velocity. 

Fully developed laminar pulsatile flow can be expected to be a superposition of a steady 
mean flow and an oscillatory flow, since the equation of motion is linear in that special 
case. For an entry length situation and/or transitional or turbulent flow, this is no longer 
true and oscillatory flow has to be treated separately from pulsatile flow. Moreover, some 
analytical solutions for pulsatile flow have singularities for zero mean velocity. Yet, results 
for pulsatile flow can possibly be taken as qualitatively representative for oscillatory flow. 
For instance the question whether the turbulence structure in an unsteady flow is different 
from that in steady flow can be discussed with knowledge about pulsatile flow. The 
majority of publ.cations deal with pulsatile flow, and these results are included in this 
survey, but listed and described distinct from those for osctilatory flow. Similarly, 
although our focus is on internal pipe flow, other geometries investigated could be a 
valuable qualitative source of information and are therefore included. From a numerical 
point of view, turbulence models for general unsteady situations may be of great 
importance and are therefore included, too. 

This literature survey extends and updates the excellent review of Seume and Simon 
(1986a) for the fluid flow problem. Since it is assumed that fluid flow is the underlying 
problem, this review does not include information about heat transfer. Once the fluid flow 
is understood, the heat transfer can be inferred. 
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1.2.2. Laminar Flow 


1 .2.2. 1 . Fully Developed Flow 

From the literature surveyed it is clear that fully developed oscillatory and pulsatile 
laminar flow are well understood. Since the Navier-Stokes equations for this special case 
are linear, a pulsating velocity can be split into a steady mean component and an oscillating 
component However, a distinction between oscillatory and pulsatile flow is necessary in 
the case of an entry length problem. 

Straight Geometries 

The oscillating flow effect was experimentally discovered for laminar oscillating flow 
conditions by Richardson and Tyler (1929), when they found the velocity distribution near 
the mouth of a pipe different than the steady state profile, i.e. the maximum velocity was 
located near the pipe walls and not in the center. Sexl (1930) was able to predict this 
behavior for sinusoidal motion in a pipe. Following are a number of analyses and 
experiments performed to study the flow pattern, pressure drop and friction factor for 
different geometries as well as for bent circular pipes: 

Analyses. Wommersley (1955) and Uchida (1956) calculated the velocities of laminar 
pulsatile flow in a straight pipe for non-sinusoida) motions of the fluid. The Uchida 
analysis is still the most prominent analysis foe laminar flow. Drake (1965) and Cedeon 
(1986) analyzed the flow panem of oscillating flows in rectangular channels of finite and 
infinite width. Drake also derived a solution for the skin friction. Vosse (1986) treated 
oscillatory flow in a flat plate channel numerically with a finite element method. Employing 
a one-dimensional analysis and the well-known steady state friction factor, Jones (1985) 
gave an analytical solution for the instantaneous pressure drop in oscillatory flow. For 
pulsatile flow , Trikha (1975) obtained the time dependent friction factor by a Laplace 
transform solution. Ohmi ei al. (1981b, c) gave instantaneous and time averaged values for 
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friction factor and pressure drop in a pipe and found behavior qualitatively similar u> 
turbulent flow. Chen and Griffin ( 1983) stated that the pressure loss in oseiltotoiy * 15 
considerably larger than in steady flow. For genera) otsseady motion. Chmbre and 
Schrock 0978) derived in inilyticll solution for fuDy developed laminar pipe How. Iguc 

„ of (1985a) obtained a dme dependent representation of the wall shear mess from the 
integral momentum and energy equations. Unsnrady bomutary layers were treated* derail 
e g by Telionis (1975). Cebeci (1986) describes an intelligent numerical scheme for 
unsteady boundary layers with large flow reversals. 

Experiments. Studies of oscillatory flow in a straight pipe ate reported by Shagal a 
al ( 1965). while Edwards and Wilkinson (1971) did the pulsatile flow case. Their results 
Show good agreement with the Uchid. analysis. Volenti 0947) perfonned experiments of 
a liquid in a U-shaped tube with free oscillations, while forced oscillations were 
investigated by Chan and Baird (1974). Oscillatory flow in tapered channels was studted 
experimentally and analytically by Cover (1986). who found that the results differed 
substantially from those in a straight channel. Similar results were presented by Ikeo and 
Uiawa 0986), who investigated the oscillatory flow pattern in an convergent tube 
numerically and experimentally. Duck (1985) solved the flow panem of a pulsatile flow m 
a nonun, form channel numerically and found that above a critical amplitude of oscillation, a 
failure of the boundary layer equations ocured. 


The flow panem of laminar pulsatile flow in bent circular pipes of various cross 
sections was studted experimentally by Chandran e, al 0 974). That of oscillating flow 
was investigated numerically and experimentally by Sudou el al (1985). Chandran el 
at. 0974) found that the maximum velocity w as shifted towards the center of cun ature 
compared to steady flow. Sudou el al. (1985) found good agreement of experiment and 
prediction A study by Sumido and Sudou (1986b) for pulsatile flow used laser-doppler 
anemometrx to measure the axial velocity in a curved pipe. They reponed good agreement 
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with their numerical predictions. Telionis (1981) gave similarity parameters for 
nondimensional treatment of pulsatile flow in curved pipes. Takami (1984) deter mined the 
pressure drop in pulsatile curved-pipe flow by a time-dependent numerical analysis and by 
an approximation method and found good agreement The wall shear stress in oscillating 
flow as a function of radius of curvature was calculated by Sumida and Sudou (1985 ). 
Yamane et al. (1985) found an increase of wall shear stress compared to a straight pipe 
under identical conditions. Finally, Sumida and Sudou (1986a) determined the pressure 
drop of an oscillating fluid in a curved pipe by experiment and analysis. 

1. 2.2.2. Entry Length Conditions 

Oscillatory Flow in Straight Geometries 

The first investigation of this kind was reported by Disselhorst and Wijngaarden 
(1980), who studied separation near the entrance of a tube under acoustic resonance. 
Thereafter separation only occurs for small and moderate Strouhal numbers. Peacock and 
Siairman (1983) predicted the entry length shorter than in steady state conditions. 
However, Seume and Simon (1986a) argued that experiments do not support this 
hypothesis. Chayrreyon (1984) proposed a time dependent entry length. Ohmi (1986a) 
measured the velocity distribution and the entry length. Apparently the only experimental 
investigation into pressure drop behavior of oscillatory flow in a straight pipe was reported 
by Taylor and Aghili ( 1984), who gave values for pressure drop which were consistently 
higher than steady state values. Their results implicitly included entry length effects. It is 
not clear how much of the reported losses were due to the oscillating flow effect alone. 

Pulsatile Flow 

Using a Laplace transform, the straight co-axial annulus was analytically solved by 
Atabek (1961), where the limiting cases (circular pipe and parallel-plate channel) received 
special consideration. 
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1.2.3. Turbulent Flow 


1. 2.3.1. Oscillating Flow 

The concept of turbulence in oscillating How is inseparable from the problem of 
transition. We will, however, distinguish between research which is primarily focussed on 
the turbulent flow structure itself, and research whose focus is the transition. The above 
mentioned experiments by Taylor and Aghxli (1983) also covered the range of turbulent 
flow. Hino et al (1983) investigated turbulent flow in a rectangular duct experimentally m 
great detail. For Re,nw*22500 and Va-309, they found that the turbulence structure of the 
oscillating flow was substantially different from steady flow, and that the accelerating and 
decelerating phases themselves had different turbulence structures. While in the accelerating 
phase turbulence was triggered near the wall but suppressed, the turbulent kinetic energy 
would be generated explosively near the wall in the decelerating phase. They also observed 
that “during the whole cycle, a layer that obeys the semilogarithmic law exists above a 
sublayer similar to the viscous sublayer”!. While in the accelerating phase, this layer was 
very thin, its thickness increased in the deceleration phase. Hino et al. also stressed the 
point that while the statistics of oscillating flow turbulence are quite different from steady 
flow, “the elementary process which maintains the turbulent production is almost the same 
as in the steady wall turbulent flow”*. In a paper in Japanese, Yoshiki et al. (1986) studied 
the velocity distributions in air between two pistons of arbitrary phase difference including 
180 c . The Re number and Va number range was 1.32X10 4 to 5.94x10* and 1 19 to 353, 
respectively. Their results indicated that the turbulence appeared “in velocity waves for all 
conditions independently of the piston phase difference. The instantaneous velocity 
distribution became almost uniform in the center region and looked like those for steady 


>M Hino, M. Kashiwayanagi. Nakayama. A and T. Ham in J. Fluid Mechan.cs, vol. 131. P 370. m 
2 m Hino, M kashi«ayanagi. Nakayama, A and T. Hara in J. Fluid Mechan.cs. vol 131. p 398, 198 


7 



turbulent flows, regardless of the phase difference” 1 . As in laminar flow, the near wall 
fluid reacted faster according to acceleration and deceleration. 

1. 2.3.2. Pulsatile Flow 

Experiments in Straieht Pipes 

Mizushina et al. (1973) measured velocity profiles and turbulence intensities and found 
that there are two distinct classes of flows according to the pulsation frequency: For lower 
pulsation frequencies, the velocities behave quasi-steady and the turbulence intensity does 
not pulsate. For higher pulsation frequencies the velocity profile is quite different from the 
steady state form and the turbulence intensity fluctuation pulsates oppositely to the velocity. 
Cousieix (1979) found, that pulsations do not significantly influence the boundary layer 
and the turbulence structure. He therefore applied a steady state turbulence model to 
simulate the boundary layer numerically. However, only a single frequency was 
investigated. Kirmse (1979) compared his own experimental data with the computations of 
Vasilev and Kvon (1971) Finding poor agreement. 

Ohmi ei al. (1980a, b,c; 1981 a,d) derived 4 characteristic parameters to describe the 
flow pattern, but without physical interpretation. They examined the influence of the 
pulsation frequency on flow pattern and turbulent friction losses. Three flow regimes 
(quasi-steady, intermediate and inertia dominant) are reported as a function of frequency. It 
is stated that the instantaneous friction factor was either equal, greater than or less than the 
steady state friction factor, but the time averaged friction factor was always greater than the 
one for steady state. Ohmi et al. (1983) investigated the eddy viscosity distribution as a 
function of pulsating frequency with a X-wire probe. The distribution was found different 
from the one in steady flow. Tu and Ramaprian (1983) studied instantaneous velocity and 
wall shear stress when the mean flow was well in the turbulent range (Rem*: 5-10 4 ). The 
frequency was varied over a wide range. Their results showed that the time mean flow was 

Yoshiki, S Tsumura. T. Endoh and N. Takama in Nihon Kikaigakkai rombunshu. vol. 52. p. 3650. 
Nov. 1986 
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affected by pulsations when the oscillation frequency approached the characteristic 
frequency of turbulence. According to their results, neither the time mean nor the ensemble- 
averaged velocity followed the universal log-law. The unsteadiness affected the turbulence 
intensity and the Reynolds-stress significantly. They noted that quasi-steady turbulence 
models neglect this effect of unsteadiness on the time mean flow. Iguchi et al. (1985b) 
studied the structure of turbulence as a function of time averaged Reynolds number, 
frequency and amplitude ratio. They presented a turbulence model including a lag in 
response time and compared it with their experimental results. Ohmi et al. (1986b) 
measured the turbulent slug and the velocity field in the inlet region of a pipe. Mao and 
Hanrany (1985:1986) measured the time variations of the wall shear stress for small 
velocity amplitude ratios and high frequencies. Their results indicate that the wall shear 
stress varies sharply over a narrow frequency range. Iguchi et al. (1986b) observed two 
types of turbulent-slug behavior, according to whether the pulsation frequency was low or 
high. In both cases the occurrence of the turbulent slug was periodic, in contrast to the 
randomness that is characteristic in steady flows. On the basis of his experiments Bulatowa 
et al. (1986 ) concluded that the turbulent fluctuations are not altered by frequency or 
channel length. He also found that a peak in the mean velocity coincided with a minimum in 

the turbulent fluctuations. 

Numerical Analyses of Straight Pipes 

Only fully developed situations have been investigated so far. Vasilev and K\ on (1971) 
used a steady state turbulence model for pulsatile flow. Their results were not confirmed by 
Kirmse (1979). Thomas (1974) used a turbulence model which utilised the mean residence 
time of a fluid particle at the wall and predicted cross-section averaged values of the shear 
stress for low frequencies Younis (1978) used the low-Re number k-£ turbulence model to 
predict the data of Lu et al. (1973). Kita et al. (1980) proposed a fluctuating (five-layer) 
eddy viscosity model to calculate the velocity distribution and the friction factor. Their 
results are confirmed by experiment. Blondeaux and Colombim (1985), using the steady 
state turbulence model of Saffman, predicted the failure of the log-law . The application of 
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this turbulence model was suggested for low frequencies only, and no conclusions about 
the genera] validity of the model were drawn. Kebede et al. (1985) used an alternative to 
the widely used concept of eddy viscosity. They replaced the Bousinesq stress-strain law 
by a set of differentia] rate equations for turbulence stresses. Surprisingly, this low-Re 
number differentia] stress model gave superior predictions than a one-equation turbulence 
model, but worse results than the low Re number k-c model. Reddy et ol. (1985) applied a 
pseudospectra] method to investigate the viscous wall region. The amplitude of die 
pulsaoons was large. The pulsation frequency was large and characteristic of the wall 
region eddies in steady turbulent flow. ‘The mean profiles of axial velocity, fluctuation 
intensity and turbulent production rate were essentially the same as in steady flow” 1 . The 
instantaneous turbulence production rate was largest at large adverse pressure gradients, 
which agrees with the findings of Hino et al. (1983) for oscillating flow. 

Other Geometries 

Flows over a flat plate were experimentally studied by Cousteix (1982; 1985) and Cook 
et al. (1985), who also did a numerical analysis. For the flat plate, Cousteix (1982) 
obtained similar findings as for the pipe (Cousteix.1979), where the turbulence structure 
near the wall was not much altered by the pulsations. Cook (1985) used an unsteady 
boundary' layer code together with a steady state turbulence model. A comparison with his 
experiments showed fair results. Experiments of Binder (1982) for large amplitude 
pulsations in a parallel- plate channel showed that the mean velocity and the streamwise 
turbulent intensity of the flow were unaffected by the large amplitude pulsations. The wall 
shear stress lead the free stream and its amplitude was less than predicted by theory. 

1.2. 3. 3. Generally Unsteady Turbulent Flows 

Gosman (1980) discussed turbulence models for the near wall region of unsteady 
flows. Given the uncertainty of whether the law of the wall holds generally, he suggested a 


*U. Reddv. J. B McLaughlin, R. J. Nunge in Fluid Eng. Trans. ASME, vol 107, n.2, p.205. June 1985 
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systematic investigation in turbulent oscillating flow in • simple geometry. He mentioned 
numerical computarions by Kouni. for fully develop* «ciltarin, Ho*. USUI, *e low 
Reynolds number turbulence model of Jams and loonier (1973). The results of those 
computations suggested that the law of the wall does not hold for numbOT “ 4 

frequencies typical for automobile engines. l,uM and Ohnd (19830 studW the influence 
of acceleration and deceletation on velocity, she* .ness, friction factor and eddy vucostty 
experimentally. On that basis they sated . limit for the applicability of the quast-steady so* 

,o an unsteady condition. In Iguchi and Ohmi ( 1983b) the authors expand the former paper 

to frictional losses in a pipe. 

1.2.4. Transitional Rows 
1 .2.4. 1 . Transition in Oscillating Row 
flow in a Pipe 

The description of transition depends strongly on the criterion used to define transition. 
This critenon is no. necessarily consistent in til publications, nor is it always sated. 

Sergeev (1966) used flow visualization to describe transition in oscillating flow and was 
the first to give an equation relating the critical Re number to the Va number. Merkli and 
Thomann ( 1975) observed transition in a resonance tube at very high frequencies. For 
these conditions they reponed a weak vonex street outside the boundary layer. A similar 
observation was made for channel flow by Sobey (198S). He also predicted these vortices 
numerically through stability considerations. Hina el al. (1976) took hot wire 
measurements of transition. Their signals showed a laminar-like phase during the 
acceleration and a mrbulent-like phase during deceleration. Crassmonn and Tama (1979) 
visualized transition by means of an electrolytic technique and reported an equation for the 
critical Re number. Ohmi el at (1982) found a large parameter range berween laminar and 
turbulent regime and quantified their findings in a transition equation. Dijkstra (198* ) 
observed transition, but did no. state a criterion. Numerical studies were done by Cayzac er 
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al. (1985). They were able to predict the lower bound of stability qualitatively, but not 
quantitatively. In his experimental work, Seume (1988) defined as criterion for transition a 
rapid increase in the measured ims velocity fluctuations. The parameter range covered in 
his experiments was representative for Stirling engines heaters and coolers. Aside from the 
always-laminar region and the transitional region he identified also an always turbulent 
region at very high Reynolds and Valensi numbers. He verified that the critical Re number 
in oscillating flow depends on the Valensi number and described two mechanisms to trigger 
turbulence: externally caused transition due to incoming fluid, and internal transition due to 
the usual boundary layer instability at higher Re numbers. All researchers agree that 
transition to turbulence can be described by a relation of the form 

Re m<max * const. ■ Va 0 - 5 

where the constant is some number around 1000. 

Other Geometries 

Park and Baird (1970) reported transition during free oscillations of a liquid in an U- 
tube. Von Kemek and Davis (1972) predicted the lower bound of stability of Stokes layers 
on a flat plate. They, like Caycac et al. (1985), could predict transition only qualitatively, 
but not quantitatively, lguchi et al. (1982) studied free oscillations in an U-tube and defined 
transition as the moment, when the velocity profile deviated from the Uchida-type laminar 
profile. Akao et al. (1986) studied transitional oscillatory flows in a rectangular duct. In 
agreement with Hi no (1976) they found that the flow had two different phases: a quasi- 
laminar one and a turbulent one. However, the flow in the laminar-like phase was quite 
different from temporary fully developed laminar flow. 

1.2.4. 2. Transition in Pulsatile Flow 

Cerrard (1971 ) probed a pulsatile flow with a mean Re-number of 3770 He found that 
closer to the wall, in the turbulent phase the velocity profile can be represented by a power 
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law according 10 (y/R)"". Sumida et al. (1984) measured pressure drop losses in 
transitional flow in curved pipes as mean-time values. Good agreement with an 
approximation is reported. Iguchi et al (1986a) studied the relanunaruation of turbulent 
slugs in a rectangular duct at four different aspect ratios. At low accelerations the behavior 
is quasi-steady. At high accelerations they observed the disappearance and reappearance of 
turbulent slugs at fixed phases in a cycle, which cannot be inferred from steady state 
behavior. Shemer et al. ( 1985a, b) found for pipe flow that the mean properties of the flow 
was not influenced by the moderate pulsations in both laminar and turbulent flow regimes. 
They presented the oscillating pan of the flow parameters as a function of amplitude and 
phase at exitation frequency. In particular, they explained the phase lag between the mean 
flow and the Reynolds stresses by the relaxation time of turbulence relative to the 
instantaneous mean shear. To capture this feature of turbulence, they proposed a complex- 
valued turbulence viscosity. Stettler and Hussain (1986) observed transition in a pipe flow- 
using laser-doppler anemometry (LDA) measurement and defined the stability-transition 
boundary as functions of Re u , mcan . Refrequency and the frequency itself. They repon 
phase-locked turbulent slugs, like Iguchi et al. (1986a). Toni and von Kemek (1986) 
examined the linear stability theory for sinusoidally pulsating pipe flow and found that the 
relevant axisymmetric disturbances are more stable in pulsatile flow than those of the mean 

flow alone. 

1.2.4. 3. Transition in Generally Unsteady Flows 

Davis ( 1 976 ) gave an extensive review for theoretical approaches to stability, which 
could be applied to oscillatory flows. Lefebre and White (1987) investigated transition to 
turbulence in a constant-acceleration pipe flow started from rest. It was reported that the 
time of transition was constant throughout the pipe, and that the critical Re number varied 
from 2x10 s to 5xl0 5 depending on the acceleration. Two characteristic parameters were 
derived to characterize the onset of transition: an acceleration parameter and a local 
boundary -layer thickness Re number. 
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1.2.5. Conclusions 


1 . ) Laminar oscillating flow is well understood. Analytical and numerical analyses agree 

well with experiments. 

2. ) There are numerous investigations on turbulent pulsatile flow, but no agreement exists 

on whether pulsations significantly disturb the steady flow pattern. Mizushina et al 
(1973), Ohmi et al. (1983) and Tu and Ramaprian (1983) reported differences 
compared to the quasi-steady pattern, while Binder (1982) and Cousteix (1979, 1982, 
1985) found no significant differences. It is believed that much of this confusion can be 
attributed to a reluctance by the various researchers to use consistent dimensionless 
similarity parameters to classify their investigations. However, the collective findings 
seem to indicate that the mean and instantaneous flow parameters are significantly 
affected by moderate to large amplitudes at high but not too high frequencies of 
pulsation. 

3. ) For transitional and turbulent oscillatory flow, especially the works of Hino et al. 

(1983) and Seume (1988) provide a pool of useful experimental data. Additional 
qualitative information can be found foremost in the papers on pulsadle flow by Tu and 
Ramaprian (1983) and Shemer (1985). So far, no numerical investigation of 
transitional and turbulent oscillatory flow in a pipe of finite length has been made. Even 
for the fully developed situation, only one investigation was mentioned ( Gosman , 
1980). 

4. ) A turbulence model which is well suited for unsteady situations has not been identified 

yet. However, it was mentioned by several authors that the turbulence model sought 
should provide a means for the relaxation time of turbulence ( Shemer et al. (1985), 
Kebebe et al. (1985), lguchi et al.(1985b)). 
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1 . 3 . Objectives 

1 . ) Develop i numerical algorithm suitable for solving the governing equations in an e^ 

and efficient manner. 

2. ) Identify a turbulence model which has the capabilities topredkt unsteady turbulent 

flow. 

3. ) Evaluate the performance of the turbulence model chosen for its capability to predict 

transition. 

4. ) If necessary, modify or optimize the chosen turbulence model. 

5. ) Compute the fluid flow and heat transfer of a number of data points representative for 

Stirling engine heaters and coolers and compare the predictions with experimental 
results if available. 

6. ) Answer the following questions: 

• Are steady state correlations appropriate representations for the fricoon coefficient and 
Nusselt number of the fully developed turbulent flow? 

• Do entrance length effects play a role? 

• Is heat transfer or fluid flow the major contributor to irreversibilities in the cases 
considered? 

1.4. Outline of This Work 

To stan with, the problem will be described in a general mathematical way. In order to 
estimate the limitations of the approach chosen it is important to state and introduce the 
assumptions made, which is done next. Then, the choice of the turbulence model is laid out 
and the model and some alternatives are discussed. At this point, the complete system of 
equations is established, whose solution should lead to the desired results. This bnngs us 
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to the numerical method for solving this system of equations, which is presented and 
discussed next. Following this, the computational results are given: First, we will test the 
numerical model and the turbulence model with known steady state cases. Then, the 
performance of the turbulence model for transition predictions is shown. Finally, the 
results of the oscillating flow computations are given in the order fluid flow, heat transfer 
and irreversibility. 


16 



2. PHYSICAL SITUATION AND ITS MATHEMATICAL 

description 


2.1. Differential Equations 

The governing differential equations express the conservation of mass, momentum and 
energy for a continuum and are given in the following for an infinitesimal control volume. 

Conservation of mass: 


+ div(pu ) = 0 

Ol 


(2.D 


Conservation of momentum: 


* 

+ pu*grad(u) * ■ grad(p) + div(T ) + f (2.2) 

Conservation of energy as enthalpy equation (neglecting radiation): 

p^l + pu*grad(h) = - div(q) + + u*grad(p) + (2.3) 

where "u is the velocity vector, p is the density, p the thermodynamic pressure and p is 
lhc dynamic viscosity, t denotes the stress tensor and * stands for any additional body- 
forces h is the specific enthalpy of the fluid and *q represents the heat flux vector. O 
denotes the viscous dissipation function and is defined as 


pd> — x*gTad(u) 


(2.4) 


It may be pointed out here that this set of differential equations together with the 
boundary conditions of Chapter 2.4 specify the problem completely, even in the case of 
turbulent flow. However, a practical solution of these differential equations involves some 
averaging ( e.g. ensemble averaging) over time and/or space. It is this process of 
averaging, during which the new unknown terms are created, for which a turbulence model 
needs to be employed later on. 

2 . 2 . Basic Assumptions 

The heat transfer and fluid flow problem in Stirling engine heat exchangers combines 
several different problems, some of which are not of primary importance of this research. 
For instance, a Stirling engine heater usually consists of a bundle of bent circular tubes. 

The fact that the tubes are bent further complicates the underlying situation but is not 
essential in order to reach the goal of this research. Therefore, the physical domain in this 
research shall be a straight tube. Also, we will assume that no changes happen in the 
azimuthal direction and that the velocity in this direction is zero, 
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Throughout this work the fluid will be treated as a Newtonian continuum. Together 
with the Stokes hypothesis, we can express the divergence of the stress tensor of eq. (2.2) 
as 


div(t ) e 2 div(ji def(u)) - -2. grad(p div(u)) 
where def( u ) denotes the deformation tensor and is defined as 
def(u) = [grad(u) + (grad(u)] T ] 


( 2 . 6 ) 


(2.7) 


4 
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Here the superscript T denotes the transpose of the tensw. With this input nnd the 
assumption of no body forces, the conservation equation for momentum becomes the 

Navier-Stokes equations: 


p + pu.grad(u) « -gradCp) - div(p grad(u)) 


♦ div[p lgrad(u)]T • \ grad(p div(u)) 


( 2 . 8 ) 


Defining the pressure as 


p = p + p div(u) 


(2.9) 


equation (2.7) can be rewritten in the form 


a— ) 1 

«£. + pu.grad(u) = -grad(P) + div(n grad(u)) + div [p Igrad(u)] J 
ch . — — 


(2.10) 


It may be noted here that, for the laminar case, a more convenient form of P could have 
been defined as 


p = p-Ipdiv(u) (2.11) 

I, can be shown that, for the case of constant dynamic viscosity, with this formulation 
the momentum equations can be expressed like eq. (2.10) but without the last term. Often, 
in laminar flow problems, the viscosity can be treated as constant. However, tn case of 
turbulent flow, the effect of turbulent mixing is frequently expressed by the concept of a 
turbulent viscosity which is not a fluid property but depends on the flow conditions and is 
therefore not a constant. In this case, the definition of P as in eq. (2.9) is preferred and eq. 
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(2.10) is applicable as written. The formulation of eq. (2.10) in axisymmetrica) coordinates 
is given in the appendix. 

Generally, the physical properties of the working fluid should be temperature 
dependent. Since we are concerned with turbulent oscillatory flows, turbulent diffusion 
will be the dominant effect, and a variation in the molecular diffusion coefficient due to 
temperature effects may be secondary. For now we will restrict the computations to 
constant properties. Once the basic problems in turbulent oscillatory flow are understood, 
variable properties may be introduced. The numerical program is perfectly capable of 
handling variable properties. As a consequence of the assumptions of constant properties 
the convective heat transfer problem is decoupled from the fluid flow problem and the laner 
can be solved first. 

In the energy equation (2.3), the first term on the left hand side may be simplified using 

Fourier's la*, 


q = -k grad(T) (2 .12) 

To transform h to T as the variable on the left hand side of eq. (2.3), we use the 
thermodynamic identity 



D J2L + 1 ( ] -pT) — 
p Dt P p 1 Dt 


(2.13) 


where D/Dt represents the substantia] derivative and 3 denotes the isobaric coefficient 
of compressibility 


|}sp 



(2.14) 
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which is zero in ease of an incompressible liquid and 1/T for u idol gas. Using the 
vector identiy 


. i, grad(Cp) 

X div[k grad(T)J * div[X grad (T)] + k — 5 — gr»d(T) 
c p p c p 


(2.15) 


the energy equation takes the foim 

p£L + pu*grad(T) = div[X grad(T)] + 7 ^ iff + u*grad(p)} * 

, Fad(Cp) 

+ k t -*— •grad(T) 

4 (2.16) 

The viscous dissipation function O for a Newtonian fluid in axisymmetric coordinates 
is given in the appendix. Assuming constant properties, the last term of eq.(2.16) is aero. 
For an ideal gas , this equation then becomes 

pg 7 pS.gradCT) - d ivljj- gradCT)) ♦ ± ■ ♦ 5«pad(p)l (J )7) 

Equations (2.1), (2.10) and (2.17) provide four equations for the four independent 
variables are u, v, T and P. They are the complete set of equations necessary to describe the 
fluid flow and heat transfer in oscillating flow conditions. For irreversibility considerations 
entropv comes into play but no additional differential equation needs to be solved. For a 
single phase single substance our thermodynamic system has two degrees of freedom. \* e 
have already specified p and T. Thus s = s(P,T), and we can solve locally for s. The 
information content of the second law of thermodynamics is implicitly contained in the 
momentum and energy equations by virtue of the definition of the stress tensor and b> 
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Fourier's law, both of which specify directions of processes. However, the differential 
equation for the entropy is useful to determine the amount of entropy generated in a given 
control volume. But since the generated entropy is a deduced quantity based upon the 
solution of the equations above, its derivation will be deferred to chapter 1 1 . 

Summarizing, the most important basic assumptions are listed in the following table: 

Table 2.1: Basic Assumptions 

• axi symmetric situation 

• constant properties 

• continuum 

• Newtonian fluid 

• Stokes hypothesis 

• no body forces 

• Fourier's law 

• No radiation heat transfer 


2.3. Dimensional Analysis 

Two different physical situations are physically equivalent if the dimensionless 
parameters characteristical for the situation are the same. In order to obtain meaningful 
dimensionless parameters, proper scales for length, time, velocity, pressure and 
temperature have to be defined. Seume and Simon (1986a) have identified the following 
scales for the oscillatory flow problem: 

length scale *scale s R = D/2 pipe radius 

time scale ( scale s ^ time for one oscillation 


velocity scale 


“scale = u m,max 


amplitude of the mean velocity 




2 

pressure scale P sc ale « p u m,max 

lemperature scale tscale * (AT)icf 

Since for now we are concerned with relative temperatures and pressures only, we will 
subtract a reference value (T* Po) from T and p. For a situation where there is a significant 
variation in the thermodynamic pressure as in real Sorting engines, it may be advantageous 
to scale the pressure with some reference pressure and not with the velocity scale 
(Recktenwald, 1989). The precise nature of the temperature scale (AT to » y« to be 
determined. For the University of Minnesota oscillatory flow experiment, a suitable 
definition would be (AT) rc f = Twall * TV Using these scales and suitable reference values 
for the thermophysical properties, the dimensionless quantities are : 


X 

* = R 

r 

r= R 

r* cot 


u 

V 



u = 

Ufn max 

V = 

U(Tl max 



R 2 

d>= 2° 

Um max 

grad - R grad 

div = R div 


P " Po 

P = 2 

P UfT1 max 

2 

F = p + ) 


T « T ' T ° 

(AT) re f 

P 

p = 

Prcf 

P 

p = 

Micf 

4# 

II 

k 

* * 

kg^f 


With these definitions the conservation equations of chapter 2.2 can be written in the 
following way. 
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Mass conservation: 


Va^H - 1 Remax div( pu)~ 0 


( 2 . 18 ) 


Conservation of momentum: 

Vap^ + y Re max pu»grad(u) « - R ^"** jraiTP) ♦ divfri gradfu)] + div[p grad(u) T ] 

( 2 . 19 ) 


Energy equation: 


Va P Ti + 2 RCmM p “' grad(T) = p7 divf 'c^ grad<T) l 


+ Va Ec |j- + | Re mw Ec u»grad(p) + EcQ 


where the dimensionless parameters are defined as: 


Re =' 


Pref u n»max ^ 


Pref 


Va 


Pref D 2 
4 Pref 


Pr 


Pref Cpref 
k*cf 


Ec = 


2 

u ^max‘ 
c Pref (^*^)rcf 


( 2 . 20 ) 


( 2 . 21 ) 


( 2 . 22 ) 


( 2 . 23 ) 


( 2 . 24 ) 
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Other, related dimensional parameters could have been be used instead, like 


Sir * 


toD „ 4 Va 

Um.max u 


( 2 . 25 ) 


Ma = 



u n*max 


Y Ry/M T 0 



(AT)ref 

To 


Ec 


(ideal gas) 


( 2 . 26 ) 


where 



( 2 . 27 ) 


Geometric similarity is maintained if UD is identical. For sinusoidal flow, the relative 
amplitude ratio, Ar. is a related and dependent similarity parameter 


Ar = 


£ u mmax 
to L 


1 D Remax 

sinusoidal flow ••• = 2 L Va 


( 2 . 28 ) 


2.4. Boundary Conditions 

The boundary conditions for these equations are: 
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2.5. Operating Range of Experiment and Points of Investigation 

Seume and Simon (1986a) have compiled die similarity parameters for a great numb er 
of Stirling engines. Based on this information, the test rig of the University of Minnesota 
was designed to cover the parameter range of most of the Sirling engines compiled. Figure 
2.1 shows the operating range of the experiment . The two offset regions correspond to 
two different pipe diameters. Within this operating range, the points chosen for numerical 



JO Y a 100 1000 


Fig. 2.1 : Data points investigated 
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3. TURBULENCE MODELING OF OSCILLATORY FLOW 

The purpose of this chapter is to document the choice of a turbulence model for further 
investigation and to discuss some of its already known deficiencies with respect to their 
impact on oscillating flow modeling. 

For now, we will restrict our attention to constant density flows. Since the Cartesian 
tensor notation dominates the turbulence literature, it shall be adopted here. 

3.1. Background 

As pointed out in Chapter 2, the basic equations to describe the fluid flow and heat 
transfer problem are valid also in the case of turbulence. However, it can be shown that in 
order to solve a problem with the mathematical model of Chapter 2, Re’/ 4 grid points in 3 
dimensions and Re 5 ' 4 time steps are necessary. An example of the required effort to fully 
simulate steady turbulent channel flow (without any turbulence model) at a Re number of 
10 000 is given by Ragallo and Moin (1984): If the smallest eddies were resolved with four 
grid points in each direction, a total of 5 x 10* grid points and 2000 time steps would be 
necessary to get to the statistical steady state. For an oscillating flow situation, a large 
number of cycles would have to be computed to filter out the periodic steady state. It would 
be impractical and excessively expensive to solve such a huge system of equations and to 
perform so many time steps. Even if a "direct" solution of the turbulence problem was 
achieved, the vast amount of generated data would have to be treated statistically, i.e. 
averaged, in order to provide manageable and meaningful information. 

Averaging the nonlinear constituent equations instead of the results eliminates some of 
the mathematical problems associated with a “direct” solution, but introduces new, 
unknown terms for which a closure has to be found by either experiment or theory. There 
are several ways by which the constituent equations can be averaged. If we are interested in 
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the coarse structure of turbulence, a subgrid model is appropriate. Here, the governing 
equations are averaged over some small volume and time interval. The flow can be 
calculated on a grid greater than this small volume and time interval. It is the effect of the 
turbulence directly in this small volume during this time interval on the flow “seen** by the 
grid which has to be provided by the subgrid model. If we are interested only in time-mean 
values, one can average the governing equations over a long time and compute the time- 
averages of the flow quantities. This is usually called Reynolds averaging. With Reynolds 
averaging, the details of the turbulence structure are lost; only the effect of the turbulence 
on the mean flow is described. Another technique for averaging is the renormalization 
group (RNG) approach. In the RNG theory, the velocity field is first transformed into a 
wave-number space by a Fourier transform. The short wave-length modes from a narrow 
wave-vector band are then averaged and their effect on the long wave-vector modes (in 
which one is usually interested) is described by a renormalized viscosity. This process is 
repeated until all scales below a certain level of wave-lengths are eliminated. RNG theory 
can generate subgrid models or Reynolds averaged models, depending on what the lowest 
allowable level of wave-lengths should be. Details of this approach can be found in Yakhot 
and Orszag (1986). 

Since the scope of this work is to give insights about practically useful quantities like 
friction coefficient and Nusselt number, the details of the turbulence structure need not be 
resolved Consequently, a subgnd model was not considered in this research. 

Once a decision .has been made which form of averaging should be used for a given 
problem, there are a number of closure options for the new’ terms created in the averaging 
process. The remainder of this chapter deals with the consequences of the closure 
assumptions made in the chosen model. 
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3.2. Phase Averaged Governing Equations 

The concept of tiroe-xveraging becomes more cumbosome in unsteady flows. 

According to Reynolds and Hussain (1972), the time independent mean flow is found by 
time-averaging the flow over some long time. If we are dealing with cycles <e.g. pulsatile 
flow), this long time must be an integer multiple of a cycle lime. A time-periodic cyclic 
quantity is obtained by "phase" averaging, i.e. ensemble avenging at fixed tunes of a 
cycle For completely unsteady flows, the only way an average can be obtained is ensemble 
averaging the results of many experiments at fixed times measured from the Stan of the 
experiment and having identical initial conditions. Since our attention is oscillatory flow, 
the time-mean of this flow is known to be aero, and long-time averaging is obsolete. The 
periodic steady state-which is the focus of this work-can be extracted from the full 
equations by phase-averaging. 

The process of phase-averaging the equations consists of three steps. 

(1 ) decompose any transport quantity into a periodic and a randomly fluctuating pan 
z = i + z ■ (2) in sen the decomposed quantities into the governing equations; and (3) 
phase-average the equations. In the following, an apostrophe will denote a fluctuating pan 
of a quantity, an overbar will indicate the phase average over this quantity. 

With this technique one is not limited to the condition that the characteristic time of the 
flow must be much larger than the Kolmogorov time scale, (v/e) 1 ' 2 . The emphasis is on the 
condition that the time behavior of the flow must have a nonrandom structure which can be 
recovered by proper averaging. Since this is the case even for high frequency oscillatory 
flows, a turbulence model may successfully be applied to this situation. 
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Continuity Equations: 





(3.U,b.c) 


The phase averaged and the fluctuating velocity fields both satisfy the continuity 
equation independently. Since the continuity equation is linear, averaging causes no new 
terms in equation (3.1b). Note that new quantities would be introduced if density 
fluctuations were to play a role here. 

Navier*Stokes Equations: 
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du 

_ l 

dt 
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where: | 

Rate of change of the 
mean velocity 

Change of mean momentum 
due to convective transport in 
mean velocity field 

Change of mean momentum 
due to convecdve transport in 
fluctuating velocity field 


Mean driving force 


(3.2) 


Diffusion of mean momentum 


As can be seen from equation (3.2), the phase averaged transport equations for the 
mean momentum contain a term involving the unknown fluctuating velocity. If no transport 


■f 
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equation it solved for the momentum transport of the fluctuant momentum field, this tern, 
must be modeled with a tutbulcnce model. However, even if mnspon equanon. are so v 
for these terns, these oanspon equations will conuin tuber, higher onJcrtcrmr winch 
originate during the averaging of those equations andwhichmus.be modekd. The po« » 
that any type of avenging of a non-linear oanspon equation like the Navw-Stokes 
equations will lead to additional unknown terms (well-known closure problem). 


Energy Equation: 



(3.3) 


where £ denotes the turbulent dissipation rate defined later. Hero, die third t ermonthe 
LHS of equation (3.3) must be modeled as well as the pressuro-diffusion term u'j 3p73x,. 
However, according to Mansour (1989). this term is negligible. 

3.3. Turbulence Models for the Phase-Averaged Equations 


It should be emphasized here that the turbulence model to predict oscillating flow 
should give correct answers for engineering applications. It should be as simple as possible 
and may very well be "custom made" for oscillating flow in a pipe and not be applicable to 
other situations. It is nfil the purpose of this research to find a generally valid turbulence 
model for any unsteady flow situation. 


The turbulence models in question can be divided into two groups, ».e. those models 
which use the eddy viscosity concept and those which directly solve an equanon for the 
term {p uj d^l- The latter models are called stress models. In a differential stress 
model, a deferential equation is derived for each component of the turbulent shear stress 
tensor, p^/. The algebraic stress model simplifies those differential equanons 
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sufficiently that an algebraic relation can be extracted. Among the models which use the 
eddy viscosity concept are the mixing length model, the k-e model and die stream function- 
vorticity model. 


3.3.1. Eddy Viscosity Concept 

The third term on the left hand side of equation (3.2) is frequently called "Reynolds 
stresses" or "turbulent stresses” because it is responsible for the enhanced turbulent cross 
stream transport of streamwise momentum and therefore works just like the diffusion of 
momentum which stems from the viscous stresses. This similarity between the laminar 
diffusion and the turbulent diffusion-like convection is the basis of a simple, old, yet very 
successful turbulence concept, the eddy viscosity concept (EVC). Noting that for now 
duj’^xj = 0, we can transform this term as 


pu 


Bu'. 

j dT 



(ty 


(3.4) 


According to Boussinesq (1S77), the Reynolds stresses can be expressed just like 
viscous stresses, but with a different diffusivity, called the eddy viscosity p t : 


- pu’.u' 
r i j 




fair du^ 

a7 + aT 
> v 



6 

u 


EVC 


(3.5) 


This equation represents the eddy- viscosity concept, in which k stands for the turbulent 
kinetic energy, defined as 



(3.6) 
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„ d 5ij is the Kronecker delta in tensor form. The last tom in eq. (3 .5) can be absorbed 
in the pressure term without loss of |enerality>. For an arbitrary scalar transport variable ». 
the EVC can be wrinen as 


— TTTT d0 

pu ;* “o^dT 


(3.7) 


where o<t, is the appropriate turbulent Prandtl number. 

Following the analogy used here, the turbulent motion is analogous to the molecular, 
the "turbulent eddy " is the analogy for the molecule and, after Prandtl, the analogy to the 
mean free path is the so-called mixing length. From kinetic theory it is known that the 
dynamic viscosity is proportional to the mean speed of a molecule times the mean free path. 
Consequently, the turbulent viscosity can be expressed as 

A A flCl 

Vi - L scale ' 1-scale 

The problem is now reduced to finding appropriate velocity and length scales and a 
proportionality constant in order to determine n t . Once *i t is known, the Reynolds stresses 
can be evaluated with (3.5) and the Navier-Stokes equations can be solved. Some 
consequences of the EVC are discussed in section 3.3.6.. 


3.3.2. Mixing Length Model 

The mixing length model, introduced by Prandtl in 1925. is still widely used today in 
industry and shall therefore be reviewed for its applicability in oscillating pipe flow. It uses 
the eddy viscosity concept of equation (3.5) and can be summarized as follows: 

. length scale = mixing length l m ; needs to be specified from empirical information 
1 Then . the modified pressure of equ (2.9) becomes P ■ p * 2/3 Im div(u ) ♦ P^j 
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velocity scale * 



( 3 . 9 ) 


• Proportionality constant *1.0 

For a 2-dimensional boundary layer type situation, the Reynolds stresses then become: 


*JT1 




du, 

5x, 


( 3 . 10 ) 


The main features of this model are: 

+ good results for simple flows 
+ simple to implement 
+ economical 

for complex flows, it is difficult to specify l m 
. does not take transport of turbulence into account 
- not suitable for rapidly developing flows and recirculating flows 

Because of the shortcomings outlined above, the mixing length model cannot be used in 
this study. 


3.3.3. Comparative Computational Tests of Various Turbulence Models 

Presently, there are many versions of k-t models and stress models available, and it is 
not obvious which model is best suited for the given task. However, most recently, a 
number of researchers conducted exhaustive tests of various turbulence models for the 
Reynolds-averaged Navier-Stokes equations. 
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P,,el Rodi and Scheurer (1985) etuunined and reared 8 nnbulenc. models for mAoleot 
flow pm. a flat plate whh and without pressure gradient. The models investigated were the 
It-e models of Chien (1982). Dutoga and Michard (1981). Hassid and Poreh (1978), 
Hoffmann (1975). Urn and Bremhors. (1981). Uunder and Sharma (1974). Reynolds 
(1976) and the k-o> model of Wilcox and Rubesin (1980). All those models are based on 
Uk eddy-viscosity concept. Patel ti al. found that the models of Chien. Lam and 
Btemhorst, Launder and Sharma and WUcox and Rubesin gave comparable results and 
were decidedly bener than the others. For tutbulence model modifications it seems 
desirable that a model bears an immediate relationship with physically measurable 
quantities This, however, is not the case for the Wilcox and Rubesin model. 

Henkes and Hoogedoom (1989) did a similar performance evaluation of turbulence 
models for natural convection flows along a vertical flat plate. They investigated the k-e 
models of Chien (1982). Hassid and Porch (1978), Hoffman (1975). Jones and Launder 
( 197 S, Lam and Bremhors. (1981). Reynolds (1976) and To and Humphrey (1975) as 
well as tire high-Reynolds number k-e model. Additionally, they investigated the algebnuc 
stress model of Cebeci and Smith (1974). Their findings were similar to those of Patel et. 
al. (1985): Overall, the models of Chien (1982). Jones and Launder (1972) and Lam and 
Bremhorst (1981 ) showed the best results. It is remarkable that the algebraic stress model 
by Cebeci and Smith gave significantly worse predictions than the k-E model above. 

Martinuzzi and Pollard (1989) compared turbulence models for steady turbulent fully 
developed pipe flow a, Re numbers of 10000. 38000. 90000 and 380000. They compared 
the high-Reynolds number k-E model, the Lam-Bremhorst low-Re number model and four 
variants of an algebraic stress model (the ones of Uunder et. a). (1975) and Naot et. al. 
(1970) both for with and without wall functions). They showed that the low-Re number k 
E model gave die best results and claimed that the use of algebraic stress models should be 
confuted to h.gh Re numbers or regions where there is only moderate shear. 
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In a 1985 paper, Kebede, Launder and Younis investigated whether a differential stress 
model applied to pulsatile flow yielded better results than the conventional EVC. They 
conclude that the EVC with k-e model has the best performance near the wall. A one- 
equation EVC model lead to rather large phase leads of uV as compared to the 
experiment, while the differential stress model produced a too large phase lag. 

Conclusions. On the basis of these tests, stress models did not yield superior results 
to models using the eddy viscosity concept Moreover, since stress models are also more 
complicated to implement and more expensive to run, a stress model shall not be used in 
this work. The models of Jones and Launder (1972) (or the updated version by Uunder 
and Sharma. 1974), Chien (1982) and Lam and Bremhorst (1981) seem to be the most 
versatile and reliable. Of those, the Lam-Bremhorst model is the only one which uses the 
true isotropic turbulent dissipation rate itself and is the easiest to implement 
computationally. Therefore we have chosen the Lam-Bremhorst model for further 

investigation. 

3.3.4. The Lam-Bremhorst Form of the k-E Model 

The k-e model uses the eddy viscosity concept , but different scales than the mixing 
length model. Whereas it has long been agreed on that >/T represents a well chosen 
velocity scale for the large scale motion 1 , many attempts have been made to conveniently 
specify a length scale. Mainly because of simple boundary conditions the use of E as a 
quasi-length scale 2 became very common. Even though the turbulent dissipation occurs at 
the smallest scales, £ is a quasi-length scale for the large scale motion. It is defined as: 


1 h is known that the turbulent kinetic energy is contained mainly in the large scale eddies Therefore \ k 
is a velocity scale for the large scale turbulent motion. 

... , ,7 -i 

kmeuc energy _ k - k . k t _£ 

2 characi time $ £, /k L L 
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will only be concerned with the isotropic put of c: 


From here on, we 


9u\ du^ 

E * v d Xj dx } (3.12) 

1 „ u* k -t model, two additional panial differential equations for the turbulent Idneoc 
energ>- and the turbulent dissipation rate are derived by manipulate the full. time, 
dependent Navier-S.okes equations. The resulting equations for k and e contain some toms 
which must be modeled. The hypotheses and assumptions going into thts closure were 
examined most recently by Mansour (1989) and shall no. be discussed here. 1. may. 
however, be pointed out here drat if a far reaching modification of the k-c model proves to 
be necessarv for oscillatory flows, the closure assumptions themselves will be up for 
discussion in his paper. Mansour pointed out that the existing models for the e equaoon 

should be improved near the wall. 

Harlow and Nakayama (1967) were the fire, to introduce a k-c model, but their model 
dtd no. predict turbulent pipe flow well. Jones and Launder (1972, 1973) proposed a 
different k-c model which gave good results for a great number of flows. Based on die 
model of lones and Launder, a number of modified low-Re number forms have been 
proposed by various researchers. All of these different forms use the same, generally 
.peed on closure assumpuons for the exact equations for k and e. The difference between 
them is how the boundary condidons are introduced and how the wall funcoons are 
formulated The particular formulation of Lam and Bremhorst (1981) offers the advan g 
tha, no additional terms are added to die k or £ equation in case of low turbulence levels. As 
mentioned above, this version wall be adopted here for further work and ts shown sn Table 



3. 1 . For a more thorough discussion of the differences between the Jones/Launder and the 
Lam/Bremhorst models, the reader is referred to the work of Schmidt and Patankar (1987). 

The so-called high-Reynolds number model (HRN) is a special case where the 
functions fy, f j and f 2 are set to 1. However, in case of fully turbulent flow where the 
turbulence Reynolds numbers Rfc and R t are high, those functions asymptotically approach 
unity. The real difference between the HRN and the low-Reynolds number model (LRN) is 
due to the boundary conditions. 

Boundary Conditions for the HRN Model. 

Neumann boundary condition dk/dr * 0 de/dr = 0 

Dirichlet boundary condition 

Neumann boundary condition dk/dx * 0 de/dx = 0 

take guidance from the law of the wall and set the near wall viscosity 
to Pi * py+/u + where u x * Cp& 25 k 0 - 5 
Neumann boundary condition for k dk/dr = 0 
Dirichlet Boundary condition for c e = u T 3 /Ky 

Boundary Conditions for the LRN Model. 

Same as above for centerline, inflow and outflow. 

wall Dirichlet and Neumann boundary condition for k k=0, dk/dr = 0 

Neumann boundary condition for t dzJ&x * 0 


centerline 

inflow 

outflow 

wall 
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3.3.5. Evaluation of the Constants in the k-E Model 


There are 5 empirical constants for the model c^, cj, C2. Ok, o £ . In order to judge 
whether the customarily used values of those constants are applicable also for oscillatory 
flows, it is useful to examine how these constants were determined for fully turbulent flow 
(f^ « fj * f 2 = 1 ). 

n The closure for G, Eq. (3.16b). applied to thin shear layers, yields 


3 

dx 

J 



(3.22) 


Substitution of (3.22) into the definition of G, Eq.(3.16a), leads to 


G = 



(3.23) 


For local equilibrium layers the generation and destruction of turbulent kinetic energy 
are in balance: 

G = E (3.24) 

Thus c cancels and 


c 




fsT 

l k J 


(3.25) 


The square root of the quantity of the right hand side was measured by Champagne. 
Harris and Corrsin (1970) to be approximately 0.3. Hence Cjj = 0.09. 
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c 2 


Batchelor and Townsend (1948) found that for grid generated turbulence at high 


Re numbers, k is inversely proportional to the distance to the grid, 
k-l/xi or k-Yl/xi 7* const. 

The transport equation for steady flow past a turbulence grid is 


(3.26) 




(3.27) 


Together with Eq. (3.26) this becomes 


e* YUj X - 2 


(3.28) 


The transpon equation for £ behind the grid is 


7riL = -c £- 
U i ax, C 2 k 


(3.29) 


Insening Eq.(3.28) into (3.29) shows that C2 = 2. Later on, cj is adjusted to 1.92. 

pp . Near the wall Eq (3.24) holds approximately and the universal law of the wall 
mav be assumed: 


♦ ln(90 

u 


(3.30) 


u = 


, ^ y ' _ /S 
y *— u ^“v p 


(3.31 a, b.c) 


where u* is the normalized velocity and y + is the dimensionless wall distance and u T is 
the fnedon velocity 
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Neglecting the convective teims, using Eq. (3.24), and substituting v t from (3.13), the 
transport equation for e becomes 




l°- c 



+ 


fc, -c 2 1e1, 

t 1 If 


(3.32) 


Equation (3.24) alone, combined with the definition for G, Eq. (3.16), yields another 
equation for c. Assuming that near the wall the shear stress is approximately equal to the 
wall shear stress and using the definition for u t from (3.31a), we can write 


2 du i 

£BU ?a^ 


(3.33) 


With this, (3.32) becomes: 


JL 

dx. 


k* 


Or r)ll 


d u./dx, 
1 * 


2^ 


8uj/5x 






(3.34) 


Using the law of the wall (3.30), the velocity gradients are evaluated as 

8 u i _ 

dx 2 K x 2 dx 2 k x 2 

From experimental data (e.g. Laufer 1954) it is known that k • 3.5 u^2 . Then Eq. 
(3.34) can be solved for c j : 


c 




V 


k 2 3.5 3 


(3.35) 


j 
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For c 2 - 1.92. tc - 0.4 and 0 £ - 1.3. C) is determined as 1.44. 

Alternatively, from Eq. (3.25). is it known thatk-u, 2 V 0J ' H* 11 - &»■ ‘ 3 ' 35) *“ 
be rewritten as 


C ‘ (3.36) 

Rodi and Scheurer (1986) argue that the assumption about the wall shear stress leading 
to Eq. (3 33) is not realistic for adverse pressure gradient Hows. 

o k «nd o7|. The two Prandtl numbers were first assumed to be close to unity 1 . 

The n " many calcularions were performed in which the consuls were systematically 
varied The values chosen are those which we believed gave the best overall agreement for 
the flows considered However, the flows considered were all steady flows. It seems 
likely that a tuning of the turbulent Prandtl numbers to unsteady flows might yield different 
results. But in turbulent pipe flows, the most dominant terms in the k and e equation are the 
production and destruedon terms. Therefore it seems unlikely that a moderate change in o k 
and ©£ has a significant effect on the overall results. 

Conclusions. The constants are determined for fully turbulent flow in the near wall 
region for simplified equilibnum situations. Even though c, is not specifically derived for 
steady flows, based on the observations of Rodi and Scheurer (1986). it appears likely that 
this constant is affected by the unsteadiness of the flow. The values for <J k and ts £ are tied 
to steady flow experiments, but it is believed that the impact of a variation of their values is 

small. 


1 W . R0 d, (1984 ) in Turbulence Models and Their Application in Hydraulics: a State of the 
Art Re* iew , 2nd ed.. lm Assoc for Hydi. Res., Delft, P-28 
2 r Hamah: and B E. Launder (1972, in J. Fluid Mech. v.52. pan.-ft. p.619 
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3.3.6. Discussion of Some Obvious Shortcomings of the Turbulence Model Chosen 

Even though the eddy viscosity concept contains many crucial assumptions, the 
measurements of Seume (1988) clearly indicate that the turbulent transport is proportional 
to the level of fluctuations and suggest that some kind of EVC seems to be appropriate for 
oscillating flows. 

A troublesome assumption of the EVC underlying the k-E model is the stress-strain 
response time. Commonly, turbulence is seen as a cascade process in which energy is 
transferred from the mean flow field to ever smaller eddies. At the end of this cascade 
process dissipation of this energy into molecular motion takes place. The amount of energy 
dissipated depends on the large scale motion, whereas the scale at which this dissipation 
occurs depends on the molecular viscosity. Clearly, this cascade process takes time in 
reality. However, the eddy viscosity concept as used here disregards this fact. With the 
EVC, a change in the large scale mean flow causes an immediate response in the turbulent 
stresses w hich are due to the action of the smaller scale motion. For strongly unsteady 
flow s it seems to be absolutely necessary to modify the EVC in order to incorporate a 
relaxation time. Such proposals have been made in the literature. Shemer et al. (1985) 
proposed that the eddy viscosity should be a complex number and mentioned earlier 
successful computations with such a model. Iguchi et al. (1985b) proposed a model for the 
axial component of the fluctuating velocity u’,ms which took the phase lag between u m and 
u'rms into account. However, this model required the experimental measurement of the 
phase lag which is impractical. 

The eddy viscosity concept hinges on the assumption of local isotropy of turbulence, 
i.e. that the turbulence structure is locally independent of direction. It is known that this 
condition is frequently not satisfied. Especially in low Re number flows, where the large 
scale and the small scale are not far apart, the assumption of local isotropy seems physically 
questionable. Despite that, turbulence models using the eddy viscosity concept have proven 
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10 be highly successful in recent years. One reason for this could be that strongly 
anisotropic turbulent eddies have little influence on the main velocity field. • theory which 
is supported by Ytkhot and Orseag (1986). Moreover, as Ytkhot and Orttag 0986) argue, 
the effects of anisotropy might be asymptotically small and may be neglected. Rodi (1984) 
asserted that “in recirculating flows where the normal stress anil shear-stress terms in the 
momentum equations are of the same order, both terms are often small compared with the 
inertial and pressure-gradient terms so that isotropy of the turbulence model is of little 

importance"' 1 . 

Rodi and Schcurcr (1986) have found thai the predictions with the k-E model become 
rather poor for a flat plate boundary layer for the case of a strong adverse pressure gradient 
and suggest a modif.canon to the model to overcome this difficulty. They showed that the 
problem of the LRN model to predict strongly adverse-pressure-gradtent flat plat boundary 
layers satisfactorily stemmed from the near wall region. Hie problem was traced to a too 
small generation rate of £ in the near wall region which leads to an oversized length scale 
and too high turbulent viscosities. However, the results of this study indicated that a 
modification of the k-E model as proposed by Hanjalic and Launder (1980) worked well for 
adverse pressure gradient situations on a flat plate. 

Even though Rodi and Scheurer’s findings uncover a serious problem in the k-E model, 
they cannot be used directly for oscillatory pipe flow. First, the convective deceleration 
over a Oat plate does not translate easily into the local deceleration experienced in pipe flow. 
In their study, the flow over the flat plate was steady at a fixed point m space. Here, the 
cascade process is statistically steady. Viewing the flow from a Lagrangian point of view, 
the decelerated fluid panicle travels through regions of steady cascade processes for which 
the EVC applies fairly well On the contrary . during a local deceleration of the flow, the 
shoncoming of the EVC will affect the predictions directly. Second, the extension of the 
Hanjalic and Launder proposal rests on the condition that the inotanonal contribution 

lW. Rodi in Turbulence Models and Their Application in Hydraulics: a State of the 

Art Review. 2nd ed.. Ini. Assoc, for Hydr. Res., Delft, p. 30 
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dum^x of the generation term G in the e equation is (i) retained 1 and (ii) multiplied by a 
higher constant than die rotational part. In our computational scheme, the imxational 
contribution is always retained and the presence of a dum/dx quantity does not give die 
desired response to a local deceleration <hWdi. Nor do we have the freedom to just replace 
dum/dx by l/u m dun»/ch. 

However, following the general ideas of Rodi and Scheurer (1986), a plausible 
modification of the k-c model would be the addition of an deceleration or acceleration term 
to the generation term of the e equation. One definition of such a term is 


Gejccc! * £ K a 

where K a is a dimensionless acceleration parameter defined as 2 


(3.37) 


M dun/*) 
*“ pu m (t) 3 dt 


(3.38) 


Then, equation (3.15) becomes 


J l J J 


(3.39) 


where D/Dt is the substantial derivative and C3 is a new constant which has to be scaled 
against experimental data. 


1 Hanjalic and Launder (1980) as well as Rodi and Scheurer (1986) used boundary layer codes which usually 
negelct terms like 


. _ 0.545 u <k*m(0 *<ji , . . 

2 Icuchi et al (1986a) defined another acceleration parameter as K * , ( ~ * TT u m ^ ) where 

pU(l) J ® 2D 

Xq; is the ume dependent quasi-steady friction factor, h remains 10 be seen whether the use of this parameter 
would vicld bcucr results. 
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Before any modified is recommended, .he results of the unmodified «*>»■»■« 
must carefully be analyzed. Only if • clear need for » modification is indicated, should 

pursued. 


3.4. Summary 

The k-e model discussed here is based on the eddy viscosity concept end uses two 
additional nanspon actions for the velocity scale end the turbuten, length 
given in the fonn of genera, unsteady transpon equations. In tts denvanon, no spectfic 
assumptions about steady flow have been made. accounts for the convecove ™pmt of 
turbulence. V,h„e the HRS vers, on rel.es on the valtdtty of the universal la. ofme wall, 
.he LRN does not. Both me HRN and me LRN models are thomughly tested and 
shown, especially the LRN version, good results for a great many flow sttuauons. In th 
denvanon of me model, no specific assumptions have been made about me stea ness 
,he fiow. Therefore mis mode, is, in principle, applicable for unsteady si.uattons and ts 

used fo; further investigation. 
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4. NUMERICAL METHOD FOR THE SOLUTION OF THE 
GOVERNING EQUATIONS 


The elliptic partial differential equations together with the boundary conditions and the 
assumptions given in Chapter 2 give a complete mathematical description of the problem. 
However, since this set of coupled differentia] equations cannot be solved analytically, 
numerical methods are needed. Any numerical method has two basic steps: First, the 
differentia] equations must be transformed to a set of algebraic equations by discretizing 
them. Second, this set of algebraic equations has to be solved. In this chapter, those two 
steps will be outlined. 

The aim of this chapter is (i) to document the discretization method used here, in 
particular the time discretization developed in this study, (ii) to describe the solution 
techniques investigated and (iii) to discuss criteria for convergence. 

4.1. Discretization Method 

4.3.1. Genera] Discretization 

Many methods to discretize differential equations have been proposed. Among the most 
prominent discretization methods are the finite difference method, the finite element method 
and collocation method. For a more detailed overview the reader is referred to Shadid 
(1989), who gives an excellent classification of the individual methods. According to 
Patankar (1988), so far no method can be claimed to be superior. The method employed in 
this work is the finite volume method which is closely related to finite difference method. 
Here, the governing equations are integrated over a small control volume. For 
completeness the standard features of this technique will be outlined. A thorough treatment 
of the finite volume discretization technique can be found in Patankar (1980) However, the 
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validity of the method shown there is extended here to 

i will be shown in detail 


variable-density time-dependent 


The first important step for an 
is to cast all differential equations into one 


convective-diffusive situations. Ibis extension < 

efficient ****«*"»** 

general fonn. Then, only one ulgonthmm 

. .. ji dependent variables. For any scalar transport viriible, the en* 6 

needed foe vtnually all ^ ^ ^ inflow into the control volume F® u5 ** **“ 

change in a control volum inflow is the sum of the 

Of generate of tins scalar wi for thc generalized scalar 0 m an 

convective and the diffusive inflow. Then, 

infinitesimal control volume 1 : 


*£^ + div(pu0-r grad<>) = S 

dt 


(4.1) 


h scalar 0 However, without loss of generality we will take the freedom to 
l ^ C . | r. u an rf «de of the equation into this source term. 

, a r anrf S which are given in Table 4.1. 
appropnate quanoues for 0, T and S w me 


if .H.c rhanie* *e will drop the overbar for the phase 

«■ L be ,0 the phase av,ta t ci ,uanut> 


averaged quantities Unless otherwise specified. 
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Table 4.J: Interpretations of 0, T and Sfor the governing equations 


Name 

Equation 

* 

r 

S 

Continuity 

2.1 

1 

0 

0 

x-Momentum 

2.10 

u 

Meff 

dP T 

-JJT* div(He!T lgrad(u)) T ) 





dp T 

r-Momentum 

2.10 

V 

Meff 

“ dT 4 divOteff Igrad(v)] ) 

Energy 

2.17 

T 

1 + Hl 

c p OT 

♦ T? grad(p) «► ♦ pc 

Turb. Kin. Energy 

3.14 

k 

Hi 

Ok 

pG-pc 




Hi 

E E 2 

Turb. Diss. Rate 

3.15 

c 

o c 

Clf]pG--C2f2p'j[' 


Now , the calculation domain fl is divided into a number of finite control volumes 
which constitute the computational grid. This grid may be different for each dependent 
variable. As shown by Patankar (1980) it is advantageous to use a "staggered grid" for 
each velocity component and a "main grid" for all other variables. The values of the 
dependent variables will be evaluated at the center of their respective control volumes. 
Figure 4.1. shows a typical grid and a typical control volume cluster and gives the 
nomenclature for the following derivations. Denoting the flux vector of the variable 0 for an 
infinitesimal control volume interface as 7 , equation (4.1) may be rewritten as 


3(pC>) 

dt 


+ div(J) = S 


(4.2) 
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v velocity control volume 


main grid point 




Figure 4.1 : The staggered grid and demils of she epical control volumes 
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This equation may be integrated in space over the finite control volume 1 and in time as 


n c t+A 


JJ J dt *** ^ / JJ div CJ)dx , «* * r’dt’* j jjs dx'idrdt' 

*W| % % w I I w 


(♦At n c 


I+Ai n i 


(4.3) 


Now we express the time integral over any quantity z as the product of the mean value 
of z during the time step times the length of the time step At, 


l+Ai 

j z(t) dt’ * z At 

t 


(4.4) 


To evaluate the individual terms in equation (4.3) a number of profile assumptions must 
be made with respect to the variation of the variables to be integrated. First, we assume that 
the mean value z is some “mix” of the old time value and the new time value, 




(4.5) 


Here, f is a "time integration factor" which may be different for each dependent variable 
and for each control volume. For now we will assume that f is a constant for each 
dependent variable. Later we will present a scheme which uses locally variable time 
integration factors. 


Second, for the first term on the left hand side and for the right hand side of equation 
(4.3) we suppose that the quantities p$ and S are constant over the control volume. Then, 
the integrated form of equation (4.3) can be written as 


* Note that the third dimension of the control volume is set io unit) in the 2 -D formulation Therefore. 

dV»dA dx where dA*rdr 
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ppQp • PpOp AV + f u + (l-fj) if * fjS + O-fj) S° <4.6) 

At 

when f , is the time integration factor for the* equation and S denies *c «*"* 
integrated scarce ten*. The superscript ° denotes the known values .. ume h • «*«“ 
absence of a superscript indicates the new rime level t*At The term D is de 

IJ m ) t ■ J w + J n * J, (47) 

JJO is defined similarly and indicates the net outflow of the integrated fluxes of the 
variable 0 time t. 

to rite same manner, the continuity equation may be in.egn.ed with mother, still 
arbitrary time integration factor f2- 

pP'ff. AV ♦ f, IF + 0 -f 2 ) IF° = 0 (4.8) 

At 1 

where 

IF e F e - F w + F n - F s (4.9) 

Here, IF and IF 0 stand for the sum of the integrated mass fluxes at the time levels 
t+At and t, respectively. 

According to Patankar (1980) we define: 

VF e 0 P = a E «V0 E ) l410) 

a £ = D e A(IPe e l) + maxl-F e ,0] (4.11) 
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(4.12) 


•! - D£ A(IP£l) -unaxI-^OJ 

J w - F w $p ■ • (0p ■ 

^ - D w AOPe* I) ♦ max[-F w ,0] 

* D° A(IPc! I) + max(-F®, 0) 

F. 

Pc, * ^ , i * e, w, n, s 

i 

A(IPc,l) = max[o, (1-O.llPe jl) 5 ] ,i s e, w, n, s 


(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 

(4.19) 


where D is the integrated diffusion flux across a control volume interface and Pe is the 
Peclet number associated with this interface. 

The definitions for a\ , *k°, as and as 0 are analogous. Furthermore it is 


^PpAV 

** “ At (4.20) 

o P> V 

~ At (4.21) 
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s e S C AV + S p 4>p AV 


( 4 . 22 ) 


For simplicity, we write 

Za NB ' a E + a s + a N + a s 


( 4 . 23 ) 


In order to transform equation 


(4.6) into an equation for the dependent variable *. 


equation (4.8) is formally multiplied by some "mix 

f 3 4> p + 0-f 3 ) Op 


of the dependent variable 0 


( 4 . 24 ) 


„d then subtracted front equation (4.6). Then, after some algebra, the disotdiabon 
equation for 0 for the control volume around P becomes 

U + fil^a^ + XF-SpAVJlOp* 

f, S<a N . B 0 OT ) * [< • d-f,) (Hb ^ S > V H *' + 
fj s c iV * (l-f,) [s°ov + £(>kb*kb)1 


(4.25) 


Note dial equation (425) does neither depend on ft nor on ft, and that no assumpnons 
were made about the values of the fs. A value of 0 corresponds to a fully expltct nme 
marching procedure (FE). a value of 1 to a fully implicit scheme (FI) and a value of 0.5 to 
the well-known Crank-Nicholson scheme (CN). Since equation (4.23) ts tndependen. o 
f2 , we can se, ft in equation (4.8) to unity without loss of generality and substitute for IF 

= - a t . For the fully implicit scheme we define 


a p.n * Sa NB + »? ' V V 


(4.26) 
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bpj-s^v+tft; 

Then equation (4.23) can be transformed to 


(4.27) 


( f i + *J *p - I P (*„,♦„,) ♦ b n ] + 

<WKie-*‘p» * (S e °*sp*^4V + «• - IF°>»°] (4 2g; 

From this is can be seen that any non-fully implicit formulation can be interpreted as a 
deviation from the fully implicit case. 

The choice of the value of fj is dictated by accuracy and stability considerations. While 
the Crank- Nicholson scheme is the most accurate, the fully implicit scheme is the most 
stable. The Crank-Nicholson scheme is mathematically unconditionally stable but can lead 
to physically unrealistic oscillations. Therefore, our goal is to develop a scheme which is as 
close as possible to a Crank-Nicholson scheme, but also gives always physically realistic 
results To get the limit of stability for the Crank-Nicholson scheme, an analytical 
perturbation analysis can only be performed for the simplified case of constant coefficients. 
In such a situation for two dimensions and an equidistant Cartesian grid, it can be shown 
(Roache, 1968) that the physical limit of stability is given by the von Neuman analysis as 


ctAt aAt 1 

Ax 2 Ay 2 < 2 

(4.29) 

uAi + vAt. < j 
Ax Ay 

(4.30) 

uAx vAv 

+ — “ < 4 
a a 

(4.31) 
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where a * I7(pc p ). It is evident that for an ever finer grid condition (4.29) is the most 
stringent condition, since the the time step must be reduced proportional to Ax* and Ay 2 . 
Computing a situation with wall turbulence requires an extremely fine grid near the wall to 
properly resolve the steep gradients there. In such a situation it is therefore practically 
impossible to apply any scheme other than the fully implicit one if a vast number of time 
steps is to be avoided. For a situation with varying coefficients and a nonuniform grid, 
another way to treat stability is the nile\ that all coefficients in equation (4.25) are to be 
positive. If this condition is violated, physically unrealistic solutions may arise. It can be 
shown that the term in the wavy brackets on the left hand side of equation (4.25) is always 
positive, regardless of the value of ft- However, it is also evident that the coefficient before 
0 pO may very well become negative for small values of f j. A remedy for this situation can 

be formulated and is shown next. 

4.3.2. Adaptive Time Integration Scheme 

As pointed out above, the task is to have a time integration scheme which at each grid 
location is as close to a Crank-Nicholson scheme as numerical stability considerations will 
allow, . The condition for stability is that the coefficient in front of 4»p° is nonneganve. Due 
to the convective-diffusive formulation of the coefficients of the neighboring points which 
we have adopted here, it is made sure that they are always nonnegative. It can be shown 
that the coefficient in front of dp is always nonnegative provided S P is formulated properly. 
The situation is more complicated for 

An equation equivalent to (4.25) can be derived in the same manner as shown above if 
one individual time integration factor is introduced for each conirol volume center and one 
for each interface (see Figure 4.2). 


] S V. Patankar in "Numerical Heat Transfer and Fluid Flow’ , Hemisphere 


Publ. Co., Washington. I9R0. 


P 
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Figure 42: Placement of time integration factors for adaptive time integration scheme 


For brevity, we will drop the index 1 in the following and define the subscript nb for 
the control volume interfaces e, w, n, s. Then equation (4.25) takes the form 

[ a i + ( ^^nb a KB + * fp^pAV } ] $p * 

I(f„ b • N ,e NB > * [>? • { 2(1 -f„b tf N D - 2(1 -urf- (1 -fp)S> v ) 1 - 

fpS^v * (l-fp)S^v * Id-U, (4.3 
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„ may be pointed ou, hem ft. inodeno — ^ace-Hux 


interface time integration factors are 


essential. This may be demonstrated in the following 


example: 


Consider two 


adjacent control volumes around the points W and P as 


shown below. 


interface flux f w J w + 0-fw) J w° 

The integrated flux across the w-interface is f* J* + 0-fw> J*° ■ 
regardless whether it is evaluated from the W or Urn P control volume. If. 
instead, only fw and fp were employed, the flux at w equated from control 
volume W would be equal to fw J. * O-fw) V while the flux evaluated 
from P would be f P l» - (Mp) Jw» • "consistency could lead to 

physically unrealistic solutions. 

The next task is .0 optimize .he individual time integration factors fp and f*. We define 
the coefficient in front of 0°. aoo. as 


. a? - I(l*f„b)a°KB 4 ^ 1 *^ )F °* °* fp)S ^ V 


(4.33) 


where f nb and f P are given an initial value of 0.5 coresponding to the Qank-Nicholson 
scheme. Subsequently, for each control volume. is evaluated and f* and f P wdl be 
corrected if a* 0 is negative. The correction sought is 


y a f a o y Af P° . Af„S^V = maxI 0 ,-a „ 

= i ^ 1 nfc a NB rbnb “ P ^ 0. m.n»! 


. miu a 1 . evaluation 
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so that when a*° is negative upon its first evaluation, it will be set to zero, and fp and 
f n b will be adjusted c o rr e spondingly. When a$° is positive, no action is taken. To evaluate 
the corrections of the time integration factors, we will assume that the correction applied 
will be the same for all time integration factors associated with one control volume, Af n b * 
Afp. Then, we can solve equation (4.34) for Af as 


Af p * Af nb 


maxfO, -ao . . . , , . ] 

6, miuaJ evajutugr 

^nb + £f - s ^ v 


(4.35) 


The AT s are evaluated for all control volume centers and interfaces in fl This will 
usually lead to a multiple evaluation of the interface time integration factors. For the actual 
correction of the time integration factors, which is of the form f„ew * finiiial + Af, the 
largesi correction will be taken. This will ensure that in no control volume the stability 
criterion is violated. Based on our experience, the CPU time increase due to this scheme is 
insignificant. 


It has been argued that the interface time integration factors can be replaced by the 
control volume center factors. Using 


£ fnb a\B = fp I a\B + I(f n b - fp) as’B 


and I f nb F « fp I F + £(f n b * fp) F 

one can asses the qualitative impact of such a measure by transforming equation ( 4 . 32 ) 


■f 
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y( a t + f F (la, |„ + EP • SrAVII 6. + H(f-v ■ V *** * ^ f nb ' V ^ jf Z 

fp^NB^NB + K • (»* f P> {Ia NB + IF ' S P AV)1 *P + 
fpSjAV + (l-fp)S®AV + (l-fp)!^^ ♦ ^ f nb ~ f pX*NB+NB * *NbW 

l^- f P> a ^' hI(f ^- f P )FOl<l>C p €4.36) 

It can be seen that the underlined terms containing the interface time integration factors 

become negligible for $p • <>£ . 0 nb “ $nb • a NB " a NB F i " ^ * <F ^ S * s 8 cnera ^ 
given for very small time steps, but for larger time steps the omission of those terms 

introduces a considerable numerical error whose consequences are not known. 

Arguing that in an unsteady situation with sufficient time steps the coefficients will vary 
little, equation (4.32) could be simplified if the old time coefficients aNB 0 wore replaced by 
the new time coefficients. While this measure seems advantageous from an economic point 
of view , this would introduce another numerical error into the discretized equations. It may 
be pointed out that the old time coefficients are readily available at the end of each time step 
calculation. Since the time integration factors also depend only on knowledge of the old 
time step, the entire term a* 0 can be efficiently evaluated at the end of each time step. 
Therefore, it is neither necessary nor desirable to employ this simplification. 

4.3.3. The Pressure Equation 

For incompressible situations, the continuity equation does not contain the the density. 
The equation of state for an incompressible fluid only gives a relationship between pressure 
and temperature Thus, an explicit equation for the pressure is still missing. In what 
follow s, we derive auxiliary equations for pressure. Furthermore, a so-called pressure 
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conection equation is derived to conect the velocity field such that it satisfies the continuity 
equation exactly. Both derivations are tied to ihe solution algorithm for pressure-velocity 
coupling. The solution algorithm employed here is based on a proposal by Patankar and 
Spalding (1972) and is of the SIMPLE type. The name "pressure correction equation" can 
be understood from the development of the original SIMPLE algorithm. In die context of 
the newer SIMPLER algorithm used here, this name is misleading since pressure is not 
corrected with this equation. For a review of other methods than die one used here, the 
interested reader should turn to Patankar (1988). 

It was shown above that the time integration scheme for the continuity equation does 
not influence the final discretization equation for <>. Therefore it may be postulated here that 
the latest velocity field shall always satisfy the continuity equation (i.e. f 2 *l). The 
discretized continuity equation is then 


a t - a? + (puA) e - (puA)* 4 (pvA) n - (pvA) t * 0 (4 y 

This equation will now be transformed to yield an equation for pressure. Independent 
of equation (4.25), (4.28) or (4.32), the discretization equation for the dependent variable 
at one point in space can be written as 


o p 0p = a E 0 £ 4 4 a N 0 N 4 a s 0 s 4 p 


(4.38) 


where the a s and P are obtained by comparison with one of the above mentioned 
equations. On this basis the u-velocity equation can be written at point e as 1 : 


o,u, - Ia nbU „ b - P + f,A,(P p -P E ) ♦ <l-f,>A,(P;.p”) 


(4.39) 


JjNoif that thf the coefficient a « represents the under-relaud coefficient if under-relaxation was done to the 
Uc equation before the pressure equauon is entered. The same applies to p and d. 
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Here the pressure' has been taken out of the source temi P and treated in the same 
fashion as all other terms contained in S c 5 . From this, an explicit equation for the velocity 

at e can be obtained 

u e = fi e ♦ f,d e (P p -P E ) ♦ 0*9 d« (PJ* p e) (4.40) 

with the definitions for the pseudo-velocity C 


fie 


E° n b U rb + g 

«e 


(4.41) 


Ae 

and the quantity de * ~ 


(4.42) 


Equation (4.40) can be used to eliminate u e and to introduce P in equation (4.37). The 
same can be done for the other velocities in equation (4.37). The resulting equation is the 

pressure equation : 


a p P p = £ (a N13 p kb^ + P 


(4.43) 


where 


<*NB = f l ( P A) nb d nb 


(4.44) 


1 Note the double meaning of "P here: the subscript P stands for the grid point, whereas otherwise P 
denotes the appropriate pressure term (cf. Chapter 2). 


2 From here on ue will restrict discussion to a spatially uniform lime 
subscript "V. A derivation for an adaptive time integration scheme is 
shown later 


discretization scheme and return 10 thr 
sifitifhi forward bui not needed as 
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®P * ^ °NB 


(4.45) 


P*a^-t,-S^ + (l- f i) (4.4l 

and is defined equivalently to equation (4.9). 

4.3.4. The Pressure Correction Equation 

After the velocity field is computed, it will satisfy the momentum equations, but not 
necessarily the continuity equation. Thus a correction to the velocity field shall be derived 
which ensures that it satisfies the mass balance exactly. Define a correction to pressure P 
and velocity u as 


P = P’ + P 


(4.47) 


u = u + u 


(4.48) 


and equivalently for v. The starred quantity denotes the quantity after solution of its 
transport equation, u' denotes the velocity correction sought and P stands for the 
corresponding pressure correction. Equation (4.39) can also be written for the starred 
velocity u*. Subtraction of this equation from equation (4.39) yields an equation for u' as 


- £ o„ b », b * f,A,(p p .p E ) * (i f,)A,(p;.p° ) 


(4.49) 


However, the old pressure field is presumed to be the known and exact; hence there is 
no correction for it. For simplicity, the first term on the right hand side is omitted. Now an 
equation for u' e can be formulated: 
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u« * fj d e ( p p * *V 


(4.50) 


Similarly, for v n ’ 

v n c f i &n ^p‘ P N^ (4.51) 

With equations (4.50) and (4.51). »’ and V can be eliminated from equation (4.4*) and 
its equivalent for v. and the resulting expressions for rhe velocities . be subsnrured mto 
the continuity equation (4.37). This yields an equation for rhe pressure correcoon P 


a p P p = X a NB P KB + P (4* 52 ) 


where 

p « a? - - IF* (4.53) 

and use is def.ned as in equation (4.43). The sum IF' is defined analogously to 
equadon (4.9) where the individual flows am evaluated with Ore starred velocrt.es. Nore 
that die right hand side of equation (4.53) is the continuity equation. In case of 
convergence, fl will tend to aero and is therefore a measure of convergence as wtll be 

pointed out later. 

4.3.5. The Velocity Correction Equations 

Based on equations (4.48) and (4.50) and knowing the pressure corrections, the 
velocity filed can be corrected according to 


u, = u; + f, d e (P p - p E ) 


(4.54) 



and 


v «* v n + ^0 , p' 


(4.55) 


4.3.6. On the Correct Choice of the Time Integration Factor for the Auxiliary Equations 

So far the pressure and pressure correction equations have been treated entirely 
equivalent to the general $ equation. This treatment lead to the factor f] in the definition 
equations (4.44) of the neighbor point coefficients and in equations (4.54) and (4.55) for 
the velocity corrections. A consequent treatment of S in the equations for u and v 
introduces fj in equations (4.39) and (4.49). It seems natural that the same values for fi 
should be used with the pressure and the pressure correction equations as with the general 
<1> equation. While implementing this scheme and testing it for laminar oscillating flow, this 
was indeed done initially. However, a thorough inspection of the predicted pressure 
showed large disagreement with the analytical pressure prediction for fully developed 
laminar oscillating flow, even though the computed velocities were very close to the 
analytical ones. The disagreement turned out to be an oscillation around the sinusoidally 
varying axial pressure distribution. A number of tests have been carried out to examine the 
influence of fj on the pressure prediction. It turns out that the best predictions are obtained 
when (i) both, the pressure and the pressure correction equation are treated as fully implicit, 
(ii) the velocity corrections are done fully implicit and (iii) the pressure source term in the 
momentum equations is treated as fully implicit (see also Fig. 4.4). This constitutes 
essentially a "staggered grid in time". As in the space grid, the velocity “time grid" differs 
from the pressure “time grid” unless a fully implicit scheme is used throughout. The code 
was implemented according to these findings. Since the auxiliary equations will finally be 
treated as fully implicit, the discussion above was formally carried out for locally uniform 
time integration factors only. 
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Ap for pipe 



Figure 4.3: Choice off] for pressure and pressure correction equation and for 
pressure source term treatment in u and v equations 
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4.2. Solution Method 


The set of algebraic discretization equations can be constituted one large matrix which 
could be solved by a direct solver. For instance Kelkar (1988) used the Yale Sparse Matrix 
Package to solve the flow and heat transfer around a square cylinder. While this has 
advantages for the coupling of the equations, this technique would be inefficient since for 
large matrices it is generally more economical to employ iterative techniques. Furthermore, 
the equations are, in general, nonlinear. Even if a direct solver was used for the entire set of 
equations, the coefficient matrix would have to be updated after each solution and solved 
again, until convergence was achieved. The intermediate solutions during this convergence 
are exact solutions only to a preliminary coefficient matrix, which is an unnecessary effort. 
Also, an iterative technique offers more freedom to treat the source terms in the equations. 
Solving turbulent flow with a k-c model involves the solution of two equations for always 
positive variables. In the course of the solution for those variables it is very important that 
their intermediate values never become negative. Intermediate negative values for k and e 
would render their solution meaningless. A technique to prevent this is outlined by 
Patankar (1980) and requires the freedom to formulate the source terms with flexibility. 
This flexibility’ is not given in a direct solution scheme. Therefore it is clear that iterative 
techniques are more suitable. 

An iterative solution method may be divided in two pans: First, the treatment of the 
nonlinearity' and the coupling technique between the individual physical equations, and 
second, the solution technique used to solve a set of linearized algebraic discretization 
equations. For the former, the SIMPLER algorithm (Patankar, 1980) was used with an 
enhancement proposed by Recktenwald (1989). For the latter, a technique was used which 
proved to be robust and most economical. A discussion of these features follows. 
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4.2.1. Solution of the Nonlinear Equations 

For the solution of the nonlinm equations. the SIMPUR elforithn. (*»*«. 1»»> 
w« used In this eljorithm. the discretieed differentW equmons ere solved sequent*^ 

With the latest velocity, density and viscosity fields, the coefficient mart, of the pressure 
equation is deteimined and solved. With the new pressure field, the new velocity field is 
computed Nest, the pressure coirection equation is solved, and, with its help, the 
velocities are corrected such that the centered velocity field satisfies exactly the continuny 
equation Ftnally. the equations for the remaining dependent variables are solved. When 
this is completed, the process is repeated a sufficient number of iterations until an overell 
convergence is reached. These iterations will be termed “nonlinear" iterations in what 

follows. 

Since the differential equations to be solved are in general nonlinear and coupled, it may 
be necessary to under-relax their solution in order to achieve convergence. Examples of 
under-relaxation techniques can be found in Patankar (1980). Kelkar (1988) and 
Recktenwald (1989). Generally, the stronger two equations art coupled, the more they 
must be under-relaxed. The degree of under-relaxation determines the speed of 
convergence of a solution. For the SIMPLER algorithm it is advantageous to have a similar 
convergence speed for each equation solved. An example of two strongly coupled 
equations are the equations for the turbulent kinetic energy and the turbulent dissipanon rate 
in the case of turbulent flow. Recktenwald (1989) observed in his calculations that when 
the solution of the pressure-velocity coupling occurred much more rapidly than the soluoon 
of the k-c coupling, the scheme diverged. It is the difference in under-relaxation factors for 
the velocity and pressure equations on the one hand, and for the k and t equation on the 
other, not their absolute level, which is responsible for this divergence. Even seemingly 
small differences in the values of the under-relaxation factors may effectively constitute 
large differences. One remedy to this problem would be to set the relaxation factors of all 
dependent variables at the lowest necessary value. This, however, would be very 
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wasteful. Alternatively, Recktenwald (1989) proposed an additional level of iterations in 
which the strongly coupled equations (e.g. k and e) are solved repeatedly one after the 
other, each time with updated coefficients. After a sufficient number of such extra 
iterations, the SIMPLER algorithm is continued. This technique was originally introduced 
especially for unsteady flow problems; it was tested during this research for steady 
turbulent pipe flow and found that in steady pipe flow this enhancement also speeds up the 
nonlinear convergence considerably. A "sufficient number" of these iterations was found to 
be typically between 1 and 5; 1 at the beginning of the computation and S near 
convergence. 

4.2.2. Solution of the Linearized Algebraic Equations 

For each grid point and dependent variable 0 an equation may be wrinen in the 
linearized form 


a p $ p « + <v> w 4 4 a s 0 s 4 p (4 56) 

The coefficients a and P can be determined by comparison with equation (4.32). For a 
fully implicit steady-flow scheme, the coefficients are identical to the ones given in 
Patankar (1980). The coefficients of the equations for all grid points for one dependent 
variable 0 constitute a matrix A, so that the problem can be written as 

A x « b (4.57) 

Many techniques for an iterative solution of this equation are available. Some have been 
tested during this work for their effectiveness. Special care was taken to ensure that the 
codes used could be fully vectorized. The codes tested were 1 

a) an unvectorized tri-diagonal matrix algorithm (TDMA) applied line by line 


'For a general dcscnpuon of the line b> line method, see e.g. Patankar (1980). 
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b) a vectorized, inverted TDMA applied line by line 

c) . vectorized, invened ted-black TOMA applied Une by Une 

d) , vectorized TOMA line by line m«hod using . OnySOUB »ba»<me » 

evaluate dot products of vectors. 

e) a vectorized SSOR algorithm 


■nte tes. problem was a sieady mrbulen, pipe flow on a 23 by 23 grid. The veetonzed 
SSOR algorithm proved .0 be extremely sensitive to the comect cho.ce of the over- 
relaxation factor. For an over-relaxation factor greater than 1.2. the aolnoon (bverg 
independently of how accum.ely.he equations wem „,v«d. The Cher four methods gave 
the following performance in CPU time based on the first method: 


8)10095: b) 45.59c c) 66.7% d) 81.2% 

•nte difference between b). c) and d) can be explained by the varying influence of the 
relative shon vector lengths (21 elements). Method b, and c) should break even when me 
vector length is greater man 128, given me present vector length of 64 words on me Cray 
2 compute,. Similarly, the use of me SCILIB subroutine pays only for very large vector 
•nte accuracy' of me solution achieved after a fixed sweep through me domarn was hr, e 

for method a, and about me same for b). c) and d). Due to me coupling and -Unease 

1 o limit Ml accuracy is needed for an intermediate solution. With 
the physical equations, only a limited accuracy is new 

this in mind, method b) was used for further work. 

Other, more sophistreared me.hods like preconditioned eonjugare gmdien, methods 
,„gh, lead to increased efficiency and accuracy. However, me optimization of me soluuon 

algorithm was not the subject of this research. 

How accurately shall the linearized equations be determined? When can me ovemll 
solution process be terminated" These questions will be discussed next. 


1 measured in terms of the residual o mlx which is defined below 
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4.2.3. Linear, Nonlinear Residuals and Convergence 

The residua] of the linearized discretization equation (4.S6) can he expressed as 



(4.58) 


where aNB k , ap*. and ^ are the coefficients af the linearized equations at a nonlinear 
iteration k, 0nb’. 0p* are the values for the linear iteration l within the solution algorithm. 
7 p is the residua] corresponding to k nonlinear and 1 linear iterations. However, the 
absolute value of r£* does not give a sure determination of whether the residual of an 
equation is small. We require that the value of Fp be small compared with otp 0^ 
(frequently the largest contributor in the sum on the right hand side of equation (4.58)). 
Therefore we scale each residual to determine its relative importance and define 


T * 1 

r P 


Jcl 


a i>0p 


(4.59) 


This scaling normalizes the coefficient matrix with respect to their diagonal elements. 
This scaling form offers the advantage, other than described in Recktenwald (1989), that 
the values at a particular point in the domain can be prescribed without rendering the 
residuals meaningless. The prescription is typically done by assigning a huge number to 
coefficient a p and assigning this huge number times the desired value to coefficient p. 

Since there is one r£ for each nodal point in Cl and dependent variable 0, we define a 
residua] vector F?p^ which contains the individual residuals as coordinates. The Euclidian 

L I 

norm R p ^ of the vector is a measure of the overall error of the computed solution of 0 in Cl, 



-kl 

K 

Po 



(4.60) 


■* 
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During consecutive linear iterations within the solution algorithm, the coefficients 
tetrnin constant and. in the limit. R " will approach aero However, it u uneconomic^ to 
drive the residual to aero at this pant because the equation solved is only the btwn«d 
form of the actual, nonlinear equation. A perfect solution for the linearised equations may 
still be far away from the solution of the nonlinear equations. The criterion to determine 
whether R * is small enough to terminate the linear iterations was 


where r£° is the so-called nonlinear residual, and 8* is some, user specified, small 
number. Recktenwald (1989) and Van Dooimaal and Raithby (1984) have discussed the 
choice of 6 0 in detail. Here, 8* was chosen between 0.1 and 0.3. 

After completion of the linear solution of one dependent variable, this process is 
repeated for the other dependent variables. 


Table 4.2: Typical Values of in the computations 


0 

u 

V 

PC 

P 

k 

c 

8* 

0.3 

0.3 

0.15 

0.15 

0.2 

0.2 


The nonlinear residual Rjj is obtained by evaluating equation (4.39) upon entering the 
solution algorithm with the latest set of coefficients. The series of the nonlinear residuals. 
R}^. from one nonlinear iteration to the next is an excellent measure of the overall, 
nonlinear convergence of a variable 0- If ^ number becomes small enough (e.g. lCri 6 ) for 
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all dependent variables, then a solution at this time step is obtained and one can proceed 
with the next time step. 

Another measure of nonlinear convergence is the maximum scaled error of die 
continuity equation, Omax* ** *°y point in Cl: 


Omax ~ 


max[abs[ a, - a® + IF* )] 

maximum flow rate across any control volume interface in Cl 


( 4 . 62 ) 


When Omax reaches a value of less than 1(H, reasonable nonlinear convergence is 
generally obtained. In fact, Omax is probably the single best overall measure of 
convergence. The nonlinear residual for the u -equation and smax are strongly correlated so 
that usually it suffices to monitor only Om*x- However, it proved useful to monitor the 
nonlinear residuals for the k and £ equations, because it is possible to reach an intermediate 
solution for which the continuity equation and the momentum equations are rather well 
satisfied but for which the k and £ equations have not yet reached convergence. 


4.3. Summary 

In this chapter, the general discretization equations were developed. A locally adaptive 
time integration scheme was developed which will be especially helpful in situations with 
wall turbulence where highly non-uniform grids are usually used. While it is more 
elaborate to implement this scheme than a fully implicit scheme, the additional CPU time 
cost is marginal. An enhanced SIMPLER algorithm for the treatment of strongly coupled 
equations was outlined. A vectorized line-by-line method was found to be a robust and 
efficient solver for the linearized equations. Criteria for the linear and nonlinear 
convergence are established and discussed. 
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PART II: ANALYSIS AND RESULTS OF THE FLUID MECHANICS PROBLEM 


s. predictions of fully developed turbulent pipe 

FLOW UNDER STEADY CONDITIONS 

In order to verify both the turbulence model itself as well as its programming, • series 
of computational tests were made to compare the predictions for fully developed pipe flow 
„ ,wo different Re numbers. One of the benchmark papers on turtulent pipe flow is that of 
Laufer (1954), where the detailed characteristics of the flow at Re numbers of 50,000 and 
500,000 are investigated Both, the turbulence model of Jones and Launder (1972) and the 
model of Lam and Bremhorst (1981) were tuned to match the set of data provided in this 
paper. To verify our solution at high Re numbets, predictions at Re = 50,000 are compared 
with the data of Uufer (1954), The paper by Kudva et al. (1972)provides data on pipe 
flow ai Re = 6000, which will be the second Re number for our test. 

5 . 1 . Predictions of the High-Reynolds Number k-e Model 

The high-Reynolds number model (HRN) computations were done with a 23 by 23 
grid with a finer mesh near the wall. The UD ratio was 150 for all computations. Figure 
5.1 show s the predictions for Re - 50 000. It can be seen that the predicted normalised 
velocity is consistently too high. But since rather few grid points were used here, this 
effect may be due in part to the grid sire. From Fig. 5.7 it is clear that the use of the HRN 
turbulence model for a Re number as low as 6000 is inappropnaie. 
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5.2. Predictions of the Low-Reynolds Number k-e Model 

The low-Reynolds number model (LRN) computations were made with a 33 by 51 grid 
with densely spaced grid points near the wall. For a Re number of 50 000, ca. 7 grid points 
were placed within the viscous sublayer (y+ S 1 1). For Re « 6000. this number was 17. 
The LTD ratio was 150 for Re = 50 000 and 250 for Re * 6000. Figure 5.1 shows that the 
computed data generally follows the law of the wall, but underpredicts u+ slightly. This is 
due to a slightly too high u T which in mm may be attributed to the relatively small number 
of grid points in the viscous sublayer. Figure 5 2 compares the computed local friction 
coefficient Cf, x for Re = 5.0x1 0 4 , and ai TI (turbulence intensity at the inflow) of 10%, with 
the experimental value for fully developed flow. Thereafter, the incoming slug flow 
develops swiftly as the rapidly decreasing Cf indicates. The minimum at about x/D=10 
corresponds to the beginning of the development of the turbulent flow structure. The fully 
developed value is slightly higher than the experimental value. Figures 5.3 and 5.4 show 
fairly close agreement of the measured and predicted turbulent kinetic energy. Figures 5.5 
and 5.6 show the comparison for the turbulent dissipation rate. In Fig. 5.6, the predictions 
are compared with Laufer's data as shown in the paper by Lam and Bremhorst (1981) as 
well as with Laufer's data as taken from the the original paper by the present author. There 
is no explanation for the disagreement of the two sources. However, the predicted data is in 
satisfactory agreement with Laufer’s data. Fig. 5.7 compares the predictions with the data 
of Kudva (1972) and the law of the wall. As can be seen , the experimental data at Re = 
6000 does not follow the logarithmic law of the wall. The measured data lie consistently 
above the log law . The predictions correctly reflect this trend. As in the case for Re = 50 
000, the predictions slightly underestimate the notmalized velocity. Finally, Fig. 5.8 
compares the predictions for the turbulent kinetic energy. Since Kudva et al. (1972) only 
report dara for u‘2, k can only be estimated from it. Fig. 5.8 shows one such estimate 
using the same ratio of u’2/ k at each radial location as in Laufer’s data. For completeness, 
also Laufer's data are shown. The prediction are in between the experimental curves. 
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Fig . 5.1: Predictions of 
fully developed turbulent pipe 

flow at Re * 50 000. 


Fig. 5.2: Computed local 
friction coefficient and steady 
state correlation for fully 
developed flow at Re-50000. 







Fit 5.3: Turbulent kinetic 
energy for fully developed 
turbulent pipe flow at Re * 
50000. 


Fig. 5.4 : Turbulent kinetic 
energy for fully developed 
turbulent pipe flow at Re = 
50000 near the wall. 




Fig. 5.5: Turbulent dissipa- 
tion rate for fully developed 
turbulent pipe flow at Re = 
50000 . 



Fig. 5 6: Turbulent dissipa- 
tion rate for fully developed 
turbulent pipe flow at Re = 
50000 near the wall. 
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Fig- 5 . 7 : Normalized velo- 
city predictions for fully 
developed turbulent pipe flow 
at Re * 6000. 


Fig. 5.8 : Turbulent kinetic 
energy for fully developed 
turbulent pipe flow at Re = 




6. TRANSITION PREDICTIONS 


In oscillatory flows at sufficiently high Remax numbers, the flow repeatedly undergoes 
transition from laminar to turbulent and vice versa. The actual Reynolds number at which 
transition occurs under those conditions is surely influenced by acceleration or deceleration 
effects and is most likely not the same as the critical Reynolds number for steady flow in 
smooth pipes. 2300. Even though our goal is to predict turbulent oscillatory flows, it is 
clear that a turbulence model which cannot predict transition for the much simpler case of 
steady flow has no promise of predicting transition in oscillatory flows. 

The objective of this chapter is to investigate the proposed turbulence model for its 
ability to predict transition in steady and in accelerated pipe flow. This chapter presents the 
results of a series of computations which examine primarily the following four questions. 

1 ) At what Re number is transition predicted for fully developed steady flow . 

2) What is the influence of the inflow boundary conditions for k and e on the 
prediction of this transition? 

3) What is the entrance region prediction for moderate Reynolds numbers (<10000) ? 

4) How do the predictions compare with experiments? 

5) How does acceleration affect the predictions? 

6.1. Experimental Observations of the Entrance Region 

It is well known that the hydrodynamical entrance length in laminar pipe flow is 

x/D = 0.05 Re 

However, there is not a consensus on the entrance length in turbulent flows. 

Commonly, we speak of entrance length as that distance from the entry along the flow 
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direction where the velocity profile changes. Wang and Tullis (1974) distinguish 3 different 
lengths: 

1 ) The length at which the wall shear stress becomes fully developed. 

2) The length at which the boundary layer reaches the centerline. 

3) The length at which the centerline velocity becomes fully developed. 

According to Nikuradse (1932), the entry length is nearly independent of the Reynolds 
number and is 40 to 50 diameters. This is in good agreement with the findings of Wang 
and Tullis (1974) who measured an x/D of 49.5. Deissler (1950) repons that the flow at the 
centerline was still developing at x/D* 100 for a rounded entrance. He also found that the 
flow close to the wall developed to its final form over a shoner distance from the entrance 
than required by flow in the center. In a later paper (1955) he quantifies this finding stating 
that the friction (factor) is approximately fully developed after 10 diameters. Bowlus and 
Brighton (1968) verify this and in addition give an analytically derived relation for the 
velocity entrance length as 


x/D = 14.25 logio (Re) * 46.0 

In one of the benchmark papers about turbulent pipe flow Laufer (1955) measured fully 
developed flow in a round pipe at a location of circa 50 diameters. And, based on his 
review of literature, Truckenbrodt (1980) states the following average relationship 

x/D = 0.6 ReO-25 . 

This overview clearly shows that the “length of the entrance region” in fully turbulent 
pipe flow is debated in literature >. It is not surprising that for transitional pipe flow even 
less is known about the entrance length of the flow. A detailed comparison of the entry 
length computations with experiments was therefore not attempted. 


J One reason ma> be that different researchers had different boundary conditions for k and e, which are not 
norma] 1> reponed or even measured 
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for Reynolds numbers up to 10000. 

6.2. Boundary Conditions at the Inflow 

For isouop.c turbulence, the turbulence intensity a. the inflow « be define. - 


T1 



(6.1) 


Umean 


The specification of a turbulence intensity a. the inflow bounds poses no problem. 
However, the spectfication of the inching turbulent dissipation nte e .s . pnob^m 
Generally, data about k and « a, inflow are no. reined b, experimenters^ Thereto*^ 
Assumptions have to be made and their impact on the predictions should be assessed. One 

option of specifying Ei n is 


vtpfr 


( 6 . 2 ) 


where X, is some function to be determined. 1. is evident that C is like » inverse 
turbulence Reynolds number Re t 


R c i = P |j E 


(6.3) 


1 See c.g Laufcr (1955). Nikuradse (1933), Deissler (1950) 
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From the definition of the isotropic turbulent viscosity we know that 


SaP k 2 


( 6 . 4 ) 


It is dear that with increasing Re also will increase. Therefore we use 

Cj, * C(Re) p 

and we define £(Re) as 1/VRe. This will give a turbulent viscosity at inflow of 


(6.5) 


p,«/Re c„p 


( 6 . 6 ) 


At a Re number of 10 6 , this will lead to a turbulent viscosity of 90 times the laminar 
viscosity, whereas for a Re of 10 000, this value is 9. Around the critical Re number of 
2300 the inflow turbjlent viscosity will be of the same order of magnitude as the laminar 
viscosity. If Re is further decreased, the numerical turbulent viscosity at the inflow will 
become less and less important compared to the laminar viscosity. This, of course, is what 
we require from physical intuition. 

Another option in specifying Cj n is to assume that at inflow, the rate of production of 
turbulent kinetic energy is in equilibrium with its dissipation rate, 


G = e 

The production rate is defined as 


du 

G 2 - u U t— ! 1 
' J dx 


p ax 

j 


(6.7) 


(6.8) 


* 
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A simple closure foT x is 



( 6 . 9 ) 


Near the wall 


du 


( 6 . 10 ) 


Defining the friction velocity as 


u T 2 = Xq/P 


( 6 . 11 ) 


we get 


G = 


Jt dx 


( 6 . 12 ) 


The model for the isotropic turbulent viscosity is 




With this, equation (6.10) becomes 


x 


. 

“ £ dx 


<6.13) 


(6.14) 


Eliminating £ with (6.7) and (6.1 1) yields 
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(6.15) 


u T -c^25 k 0.5 

Combining (6.4) with (6.10) and (6.14) yields 


£ 


: k 2 *\ 


Using the universal law of the wall, we get 


(6.15) 


cjj”k IJ 

C ’ * >'in ( 6 . 16 ) 

The choice for yj n is either one representative length scale for the inflow ( e.g. the 
radius of the tube) or the distance to the wall. This, however, will lead to a singularity of £ 
at the wall and a rather low value at the centerline. As a consequence, the turbulent 
viscosit> at inflow will peak sharply in the center and be very low near the wall. 

A third w ay to specify £ is to infer the nature of from experimental data for fully 
developed pipe flow . Thereafter, the friction velocity can be expressed as 


v f 0.197 Re 0 873 Re <41 0 4 

B _o I 

D l 0.151 Re 09 Re>410 4 (6.17) 

From Schlichting (1980), v t>max /(u-t R) * 0.09 * c^. 

With (6.17) this can be rewritten as 

r 0.1 Re 0 875 Re<410 4 

l 0.076Re° 9 Re>4-10 4 (6.18) 
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In order to have ^asa function of the turbulence intensity at 


the inflow, we define 


IV 


c u p T1 Re 


.0.175 


( 6 . 19 ) 


which gives a function £ when put into (6.2) of 


C(Re) 


T1 Re 0 875 


( 6 . 20 ) 


It is this formulation that we have adopted for the LRN computations. 

6.3. Transition Predictions of Quasi-Steady Flow 

The high Reynolds number version of the k-e model c.nnot .ecoum for any transitional 
effects. It incorporates only the turbulent viscosity and neglects the influence of the laminar 
viscosity entirely This is the major reason why a high-Reynolds number tutbulence model 

is not used in this study. 

Contrary to the high-Reynolds number turbulence model, the LRN version can 
potentially predict transition to turbulence and relaminaiirarion since it takes the effect of the 
molecular viscosity into account where necessary. It has been shown (Jones and Uunder 
(1972). Launder and Spaldtng (1974). Schmid, and Patankar (1987)) that low Reynolds 
number models ore capable of predicting transition, at least qualitatively. Schmidt and 
Patankar (1987) investigated the prediction performance of the models of Jones and 
Launder and of Lam and Bremhorst in steady flow over a flai plate. They found that the 
starting location of transition was predicted too early and that the length over which 
transition occurred was undeipredicted. Jones and Uunder (1972) rtpon the performance 
of their turbulence model for fully developed pipe flow. It can be seen that their model 
predtets transition at a too low Reynolds number of 1600 as well as a too narrow range of 
Rc over which the transitional state of the flow prevails (Figure 6.1). 
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Fig. 6.1 : Transition prediction for fully developed pipe flow with the model of Jones and 
Launder (1972 ). 



Fig. 6.2 : Transition prediction in fully developed pipe flow with the mode! of Lain and 
Bremhors t ( 19SI ). 
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k this study, two series of transition tests wera perfonned, _ f« «* . ■“ 
My. »0">«r one for 0,5*. Most of the tests wera p«rf«*d ”***“'“ 

grid points, whera the «U grid .ines w« - *• *"“* " 

densely spaced near the wall. This grid was sufficient for the nwtoate **y*M*« 
under investigation and still Prided reenable convergence raras ( as «">!»<* to targer 

grids). 

k a. following, we define transition of the fully developed flow as the point when the 

computed friction coefficient a, the downsueam end of the pipe stmsto deviate from the 

corresponding laminar value According to this definition, bod. tests gtve the same rasu . 
independent of the turbulence intensity a. the inflow and are in line with the findmgs o 
Schmid, and Patankar (1987) for oansition predictions for flow over a flat plate and Jones 
and Launder (1972) for fully developed pipe flows. Figure 6.2 shows the results of the 
friction coefficient computations as well as the measurements of Nukumdse (1932). As can 
be seen, oansition is predicted a, a Re number of 3450 which comesponds to the upper end 
of the Re-number transition range shown by Nikuradse. 

The Re number range over which transition is predicted is much smaller than measured 
for example by Nikuradse (1932. ca. 2 Recr). At one point in the sequence of 
computations, the predicted c, values jump suddenly from the laminar values of 16/Re to a 
turbulent value when the Re number is incraased slighdy. This does no. properly reflect the 
real intermittent transition process which occurs over, rather broad band of Re numbera. 

One must bear in mind that the results shown in Fig. 6.2 are obtained for UD ratiosof 
up to 500. In a Stirling engine beat exchanger, and in the experimental test rig for 
oscillating flow research a. the University of Minnesota, the UD ratio ts much less. It ts 
therefore of pea. practical interest to examine the predictions of the model with regard .o 
the developing flow. For low turbulence intensity levels a. .he inflow, the usual lanunar 
flow behavior with laminar entrance length and parabolic fully developed profile is 
predicted up to a Reynolds number of 3450. A, high Re. the flow will fust develop as in a 
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laminar flow. At some point downstream, die flow will become unstable and undergo 
transition from laminar to turbulent. Since this transition takes place over a very sm al l axial 
distance, we will speak of a transition front. One example of such a transition front is 
shown in Fig. 6.3. Further, it is interesting to note that the predicted transition occurs 
simultaneously over the cross section. The location of this transition front depends largely 
on two factors: (i) the level of turbulence intensity at die inflow (TI) and (ii) the Reynolds 
number. 

Turbulence intensity dependence. As the turbulence intensity is decreased, the 
transition from moves downstream . In the limit of a very low turbulence intensity the 
transition front approaches a maximum downstream location and does not move any further 
(see Fig. 6.4). This can be explained from a physical or numerical point of view: In 
addition to the imposed turbulence intensity at the entrance of the tube, in real pipe flow 
there are always disturbances downstream. If the Re number is sufficiently high, these 
disturbances will cause the flow to become unstable even if extreme caution is exercised to 
have a very low level of turbulence intensity at the inflow. Numerical approximations also 
act like physical disturbances. Even if we turn off the imposed turbulence intensity, there 
will be a residual level of “numerical disturbances" in the domain which will cause 
transition. A lower turbulence intensity also makes the transition front steeper, more 
abrupt. 

As mentioned above, the results for the fully developed flow are independent of the 
turbulence intensity at the inflow. However, in the developing region in the transitional Re 
number range the turbulence intensity level has a decided effect on the flow for most of the 
length. At Re « 3450, the flow well downstream will be predicted ultimately to be laminar. 
For a Tl of 0.5%, the flow follows a normal laminar behavior throughout. Increasing the 
TI at this Re number creates a region where the flow looks very much like a turbulent flow 
over much of the tube length. For example, a TI of 10% will lead to a seemingly fully 
developed turbulent profile at an x/D beyond 100, but at x/D • 250, the flow suddenly 
relaminarizes (Fig 6.5 and Fig. 6.6). This implies that, contrary to the fully developed 


j 
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Fig. 6.3: Transition from example for a Re number of 6000 and a turbulence intensity at 
inflow of 0.5 5c. 
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— i 1 
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Fig. 6.4: Influence of the turbulence intensity at irrflow on the location of the transition 
front at Re = 6000: normalized centerline velocity vs. axial distance. 
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Fig. 6.6: Entrance length relaminarization at Re = 3450 for different levels of turbulence 
i intensity at inflow: normalized centerline velocity vs. axial distance. 
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case, the developing region of the flow at moderate Re is strongly affected by the boundary 
condition at inflow. 

Reynolds number dependence. As the Reynolds number is increased, the 
transition front moves upstream. Fig. 6.7. At higher Re numbers <e.g. >5000) and no* mo 
low turbulence intensities (e.g. >2%) the transition front can hardly be seen any more. 

Then, transition takes place practically instantaneously at the entrance of the tube. This is 
probably the reason why the existence of such a transition front is not mentioned in any of 
the reviewed experimental papers on turbulent pipe flow. 

An interesting situation occurs around the nominal value for transition for the fully 
developed flow . When the Re number is decreased further, the transition front moves 
downstream. With very low turbulence intensity at inflow, the location of the transition 
front stands at an x/D of around 300 for Re « 3470. When the Re number is lowered 
further to 3460. the transition front hardly moves downstream any more, and immediately 
after the front the flow starts to relaminarize, Fig. 6.8. The final profile at this Re number 
looks very much like a normal laminar profile. Here, the transition front looks like a cut 
into a laminar flow. Low ering the Re number even further just reduces this "cut" until it 
completely disappears at Re = 3450. 

Initial guess dependence. At very low inflow turbulence intensities, the predicted 
location of the transition front is also affected by the choice of the initial guess for k and £ 
in the computational domain. Here, the initial guess for k is under investigation. A typical 
initial guess for k and e inside of the computational domain is just like the specification of k 
and £ at the inflow . For the results shown in Fig. 6.9. the initial choice of the turbulence 
intensity inside of the domain was varied while the T1 at inflow was held constant The 
results can also be explained with the action of "numerical disturbances as above. Each 
initial guess also represents an initial error and leads to a particular level of numerical 
disturbance In the near vicinity of transition, both the laminar and turbulent solutions for 
the equations are permissible and equally likely. In reality, such a situation would be 
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Fig. 6.7: Influence of ihe Reynolds number on the location of the transition front for 
059c turbulence intensity at inflow: normalized centerline velocity vs. axial distance. 
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Fig. 6.9 : Initial guess 
dependence of the prediction 
of the transition front for low 
turbulence intensity at inflow 
(0J%) and Re * 6000: 
normalized centerline velocity 
vs. axial distance. 



Fig. 6.10: Initial guess 
dependence of the prediction 
of the transition front for 
higher turbulence intensity at 
inflow (29c) and Re = 6000: 
normalized centerline velocity 
vs. axial distance. 
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Observed as imennittem. In this test we simulate a awdy state situation, and the pmgmm 
finally must decide on one particular solution. Since both solutions are equally likely, the 
decision to converge in one way or the other may be influenced by small disturbances, 
which may vety well be the ones resulting from the initial disturbance. At higher urtukeee 
intensities, this effect is practically unnoticeable (see Ftg. 610). It should be noted here 
that during the course of this research the convergence speed of the program was increased 
dramatically by various measures. This, however, changed the inherent numerical 
disturbances such that earlier results with slower convergence could not be recreated 
exactly with the enhanced version of the program. 

Domain and grid influence. In this work it was also verified that the existence and 
location of the transition front is not a grid or outflow.bound»y effect. It was found that a 
certain minimum number of redial grid points is needed to predict the front, but beyond this 
number, the front was predicted consistently at about the same location (Fig 6.1 1 ). 

Discussion. It is worthwhile noting that the transition front can exist at locations 
where the flow should be fully developed. Even though there is disagreement over the 
length of the turbulent entrance, it is clear that spatial transition at low n occurs drier this 
"entrance region". For the nominal, fully developed transition Re number predicted, the 
laminar transition length is 172.5; the transition front for Re = 3460 develops only at 
around x„/D - 300. These findings suggest that at low Re numbersand low T1 the entrance 

length is much longer than previously assumed. 

•me results of this study suggest that the transition process in a finite pipe must be 
described by two parameters. Re and Xff /D, where Xtr /D is the transition front cross 
section. Or, alternatively, the number of parameters can be reduced by adopung the notion 
of an external boundary layer and working with a momentum thickness Reynolds number. 
However, speaking of the transition is misleading because there is spatial transition from a 
Lagrange point of view, and there is ordinary transition from an Eulerian point of view 
very far dow n the pipe which depends on the Re number only. 
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Fig. 6.11: Domain and grid influence on the predicted location of the transition front: 
normalized centerline velocity vs axial distance. Top: Influence of radial grid at Re = 
3500 . Bottom: Influence of axial grid at Re - 6000. TI = 0.5%. 
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Conclusions. For developing snd fully developed pipe flow. ,u.K»dvely mi 
quantiudvely c«rec. trensition predictitms canbemwieby Ae UnvBremhres, ton. of Ae 
LRN k . £ turbulence model, but Ae mnsidond Re number mnfe to »o wnow. Even 
though the prediction to the fully developed flow is insensitive to the turbulence intensity 
the inflow, the developing region to very sensitive toil with reg*d to trensmon. Stnce 

the practical applications-which we are ultimately taerested in-will limn the snuaooos to 
mostly developing flows, we can say that for this model, transition in the developutg 
region can be triggered by the choice of boundary conditions to It and e. This is very 
desirable in light of the findings of Seume (1988) who concluded that transition in 
oscillating flow is often determined by Ac state of the fluid before flow reversal. Hus flutd 
might have been outride of Ae computational domain at flow reversal and might be entenng 

the domain during the computation. 


6.4. On the Reproducibility of Transitional Steady Flow Results 


Most of the results show n in chapter 6.2 were produced with a 33*51 grid, the same 
grid as used for Ae first computation of oscillating flow with a LRN turbulence model. It 
was shown Aat the initial guess for k and c docs have an effect on the prediction of Ae 
location of Ae transition from (Fig. 6.9) It had been found Aat a minimum number of 
radial grid points is necessary to predict spatial transition and Aat, generally, more rad.al 
grid points shift the transition front to lower a/D’s (Figure 6.1 la). However. Ae eatrtence 
of a transition front seems to be independent of grid and domain effects, as long as a 
sufficient number of grid points are taken. A close look at Figure 6.1 1 shows Aat these 
conclusions were drawn a, Re numbers of 3500 and 6000. boA above Ae predicted 

"critical” Rc number of 3450. 


Matters are more complicated ai the critical Re number of 3450: The findings of Figure 
6.5 are quite reproducible for the exactly the same conditions (i.e. initial guess, grid. T1 
etc ). However, it was found that , for a TJ of 25 and a finer grid of 63 by 63 gnd points. 
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the relaminarization of Figure 6.5b could not be reproduced. Rather, the flow would 
converge to a fully turbulent situation (Figure 6.12)! To test the stability of the coarse grid 
solution (33 by 51, Fig. 6.5b), it was transposed on a 63 by 63 grid and served as the 
initial guess for a continuation of this solution with die fine grid. The results of this 
continuation are identical to die initial guess and die solution does not change (Figure 
6.1 3). Vice versa, the results originally obtained from the fine grid solution (63 by 63) 
were transformed to a coarser grid (33 by 51) and used as initial guess for a continuation of 
this solution on the coarser grid. Similarly, the results were identical to their initial guess, 
and no spatial relaminarization was predicted (Figure 6.14). This leads to the conclusion 
that a continuation of a converged solution obtained with one grid will not change these 
converged results independently of the grid used in the continuation. Additional 
computations (no continuations) done with a 51*51 grid and a 33 by 63 grid did not predict 
the relaminarization of Figure 6.5b (Figure 6.15). 


Table 6.1: Grid influence at "critical" Reynolds number of 3450 


Re 

n 

grid 

continuation of 
other grid 

relaminarization 

predicted? 

Figure 

3450 

27c 

33*51 

- 

yes 

6.5b 

3450 

27c 

33*51 

63*63 

no 

6.14 

3450 

27c 

51*51 

- 

no 

6.15 

3450 

27c 

33*63 

- 

no 

6.15 

3450 

27c 

63*63 

- 

no 

6.12 

3450 

27c 

63*63 

33*51 

yes 

6.13 


For a T1 of 0.5^, the 33 by 51 and 51 by 51 grids predict fully laminar flow, whereas 
the 33 by 63 and 64 by 64 grid predict fully developed turbulent flow (Figure 6.16). These 
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, , t>* *A*n nnd TI s 29c Top Relaminarization predicted wiih 

fig 6J2: Predictions for Re = 3450 Qna l J r 

33 bi 51 grid Bottom: Fully turbulent flow prediction with 63 by 63 grid 
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Fig. 6.J3: Continuation of results of Fig. 6. 12a with 63 by 63 grid. 



Fig. 6.14: Continuation of results of Fig. 6. 12b with 33 by 51 grid. 












findi „,s eta* show tha, there is a severe grid influence on the Colons the "crt^T 
Re number of 3450. 

Figure 6.17 shows enoiher influence: both velocity profiles Be computed **«•* 
euctly same conditions, except for the under-reUxation factors for k «>d c. is 
that the loci of the transition front re not identical. 

The next consideration was whether qualitatively suntlr results as described ioeita® 1 
6 2 (e g at a somewhat different Re number) are ptedteted for a finer gnd f y • 

number of conations detent that the new Recr*-*: * ™ * 

2985 (Figure 6.18). However, a. Re * 2990, the flow underwent a normal spaoal 

transition a. around ./!««., whereas a, Re-2885 a oansition did no. occur. A 
relaminarization of the flow as shown in Figure 6.8 for 3460 was no. deteoed wtth th» 


grid. 

Table 62: Predictions at around Re - 2985 
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points). Two computations at Re numbers of 2985 and 2995 did not predict transition for a 
TI of 0.5% (Figure 6.19). 

At the new “critical** Re number, the influence of Tl was studied. While for a Tl of 
0.5% the flow is laminar throughout, a Tl of 10% causes die flow to undergo transition to 
turbulent without spatial relaminarization. This means that for the 64 by 64 grid the TI 
boundary condition does influence predictions for die fully developed regime. This is a 
dear contradiction to the statement made earlier, ie. that the Tl does not influence the fully 
developed regime (Figure 6.20). 

It is interesting to see how a predicted solution changes if the TI is varied suddenly. For 
this, the converged solution of Re * 2985, Tl « 0.5% (64 by 64 grid) was taken as initial 
guess for a computation with TI set to 10%. The result of this computation can be seen in 
Figures 6.21 and 6.22 as a function of the number of iterations. As one sees, the laminar 
portion of the flow is slowly pushed out as the iterations proceed, or-altematively--with 
time, the flow in the pipe will become turbulent when all of a sudden the low turbulence 
intensity at the inflow is replaced by a high one. Alternatively, if the converged solution for 
Re = 2985, Tl = 10% is taken as an initial guess for a computation where TI = 0.5%, the 
greatest pan of the flow remains turbulent. The laminar ponion of the flow near the 
entrance of the pipe is extended only slightly due to the sudden decrease of the disturbance 
level (Figures 6.21 and 6.23). 

Conclusions: 

1 ) In the transitional Re number range the results of chapter 6 are 
* qualitatively valid. 

- quantitatively valid only for the very conditions for which they were established (e.g. 
a 33 by 51 grid, specific under-relaxation factors, etc.). 

2) Transition in the fully developed regime as well as spatial transition is influenced by 
numerical disturbances which in turn depend on grid, under-relaxation, initial guess for 


u« various variables, etc.. Therefore there is no one «*W I”' 0 ™ 1 " rf 

the k-E turbulence model. 

3) For tower and higher Re numben, toe Pterions « v^id »d 

4) The performance of toe turbulence model reflects toe physics insofar >* two acltmons, 

umely the laminar and the turbulent one, become equally likely around nan 

“compete” for dominance. . 

5, For about die first 100 diameters downstream from die pipe enoance die contradicoon 

with regard of toe sensitivity of fully developed flow toll does not have severe 

consequences. This region is definitively affected by toe Has shown by -1 

computations. 

t.S. Transition Predictions of Constant Acceleration Pipe Flow 

The nex, logical step to check toe performance of toe k-c model regarding oansidon to 
turbulence is to apply i. to this physical situation: A fluid in a tong pipe is initially a. resc A, 
tone - 0, toe fluid is accelerated a, a constant rate. As toe Re number leases, toe flow 
wffl undergo transition a, a Re number higher than in quasi-steady flow. Thts MW — 
experimentally investigated by Lefebre and White (1987). Comparison of toe numerical 
predictions with toe experiment allows the turbulence model to be tested for acce eranon m 
a transitional situation. Although this experiment has much in common with toe accelemmg 
phase in oscilladng flow, there is one major difference: Here, the flow before ™ sino "‘ s 
absolutely undisturbed, and the disturbances develop only a. transition. In osctUarin* flow, 
there is some level of residual turbulence before oansition which will most probably trigger 
transition earlier than seen here. While we may expea to see some grid dependence for . 
Situadon considered in this chapter, this dependence will be significantly smaller in 
oscillating flow because of the residual turbulence. 

The goal of this study is to determine how well toe turbulence model predicts oansinon 
to turbulence in accelerated pipe flow . 
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Fig. 6 .17: Predictions for Re = 3450 and Tl = 03% for various under-relaxation 
factors a Top: a* = etc - 0.5; bottom: ai-a c ~ 0.7. 
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Fig. 6.18: Influence of the 
Reynolds number on the 
location of the transition front 
for 0.5% Tl at inflow: 
normalized centerline velocity 
vs. axial distance. 



Fig.6.19: Grid dependence 
of the critical Re number 
obtained on the 64 by 64 grid. 
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Fig .6.2 la: Convergence of 
solution when Tl is suddenly 
varied, a) Start from the 
solution for Tl * 0J%. 
continuation with Tl * l(fk. 


Fig. 6.2 lb: Convergence of 
solution when Tl is suddenly 
varied. b) Start from solution 
for Tl - lOft, continuation 
with Tl = 0_59c 
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Fig. 6.22: Convergence of solution when Tl is suddenly increased. Start of the calculations from the solution for 17*03%, 
continuation with TI=IO.O%. Numbers indicate the number of nonlinear iterations, i.e. the computational time. 






Lefebre and White (1987) used a 30 m long test section (ID * 5 cm) with water as the 
working fluid. The flow was accelerated from rest to 9 or 11 with accelerations 
between 1.85 and 1 1.8 nVs 2 . Their LDV and surface shear stress sensor measurements 
showed that the flow underwent transition throughout the lest section at practically die same 
time. An acceleration parameter was defined as 


P dum(t) 
K *' pum0) 3 * 


( 6 . 21 ) 


as well as a dimensionless time 


r* 



( 6 . 22 ) 


In their experiments, the transitions times i reported are from 0.00105 to 0.0032. The 
acceleration parameter at transition was nearly constant for the experiments at K*,u * 
1.53xl0‘ 8 . However, during the acceleration, this parameter changes due to its definition. 


Computations have been made for the same non-dimensional situation as in the 
experiment. To translate the dimensionally given acceleration data of Lefebre and White, a 
non-dimensional acceleration was defined as 


. 2 P 2 R3 dum(t) 

K » * dt 


(6.23) 


Since dU/dt was constant for each experiment, this parameter remained constant during 
acceleration. Using the propen)’ data of water and the ID of the experiment, the non- 
dimensional acceleration could be determined for different cases 1, D and ID. Case 1 is set 
arbitrarily to a very low acceleration, case II corresponds to the lowest experimentally 
investigated level and case III to the highest investigated level of acceleration. As will be 
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Table 6.3: Results of accelerated pipe flow computations 








1 

n 

ID 

experiment 

Ka* 

5.0e+6 

5.724e+7 

3.65 le+8 

same 

K« «r 

2.5e-6 

3.58e*6 

3.91e-6 

1.53e-8 

Ren rr 

2.0e+4 

4.0e+4 

7.2e+4 




7.0e-4 1 

2.0e-4 


699 

R ffr° 5 

2.37 

3.17 

3.11 

2 . 80 1 

Rcg.tr 

1500 

1680 

1584 

23200 ± 15* 


seen later, the range of non-dimensional accelerations corresponds to the range found in 
the different test cases for oscillating flow. The results of the computations are listed in 
Table 6.3. In the computations, transition to turbulence was defined to happen at that time 
step when the computed turbulent viscosity at any point in the domain became at least of the 
order of the molecular viscosity The findings of this study are. 

1. ) The turbulence model displays qualitatively the correct behavior: As the acceleraoon 

increases, Reo,u increases. 

2. ) Quantitatively, transition to turbulence is predicted about one order of magnitude too 

early ( at too high values of the acceleration parameter and at too low times 
corresponding to too low Rep)- 

3. ) The computational value of K a>u is varies only slightly for the different experiments; in 

the experiment this value is also almost constant. 

4. ) The computed non-dimensional boundary layer thickness is close to the theoretical 

value of 2.85 cited in Lefebre and White. 

5. ) The predicted transitional Re number based on the 99* boundary layer thickness is 

about one order of magnitude too low. 


1 Conclauon equation (4) given by Levebrt and White (1987; 
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Conclusions. While the turbulence model displays qualitatively the correct results, it 
predicts the actual onset of transition about one order of magnitude too early. Given the fan 
that (i) because of expected grid dependence, only qualitative results were pursued and (ii) 
the transition criterion applied is somewhat arbitrary, the performance of the turbulence 
model is viewed as being satisfactory. However, this study is important for suggesting 
possible future improvement in the chosen turbulence model. 
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7. PREDICTIONS of oscillatory flow 

Since turbulence models we generally highly empirical, their predictions should be 
compared with experiments when they are applied to a new stuadon. A navel situation » 
given when the k-e model is applied to oscillating now. and comparison with the data of 
th, oscillating now test facility of the Univemty of Minnesota will be owde. But accctttag 
» which criteria will the agreement between the computed prediction and the experiment be 
judged ? Ideally, a direct comparison of the turbulent shear stresses seems to be desirable. 
However, shear stresses have not been measured and might not be measured at all. 

Experimental data documenting the transition in oscillatory flow have been established 
by Seume (1988). This data can be taken to check the muisition predictions of the 
turbulence models qualitatively. A quantitative comparison is not possible, since (i) the 
transition data itself art valid only qualitatively (cf. Seume. 1988) and <ii> the only 
measured fluctuation component is u' . and a comparison with the quantity k of the 
turbulence model involves knowledge of how isotropic the turbulence is. From the data of 
Laufer (1954) it is evident that k# 1.5 u^ even for steady flow. 

Measured velocity profiles are available for one data point (termed SPRE in this work). 
Quandtative agreement of the predictions with the experimental data can be checked. For 
engine design considerations it is important to know the local and average fricnon 
coefficient, c ( . which is proportional to du/8r at the wall. Exact prediction of the velocity 
gradient near the wall is difficult because of the steep gradients there. Therefore, a 
comparison between the computed c ( and the measured is desirable. While in principle the 
experimental data reveals this information, due to experimental problems it has not been 
possible so far to measure cj. Computationally, C( data is available whenever velocity 

profiles are computed. 

Thus, the availability of experimental data allow the following comparisons to-date: 
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1) qualitative prediction of transition 

2) quantitative agreement of velocity profiles not too close near die wall. 

Nomenclature: In the following, the nomenclature and comparisons will refer to the 
University of Minnesota oscillating flow experiment. For a description of the experimental 
set-up see Seume (1988). The the x-axis starts at die “drive end" of die tube and ends at the 
“open end". Four out of the five cases investigated here used a pipe length L/D of 60. 
Consequently, an axial location x/D of 44 is closer to the open end and an x/D of 16 is 
closer to the drive end. The mean flow velocity will be of the type 

Um(t) * Um .max sin(CDt) (7.1) 

The cycle time will be expressed in terms of crank angle. However, contrary to Seume 
(1988), here we define crank angles between 0° and 180° when the mean flow is along the 
positive x-axis (i.e., coming from the drive end), while crank angles between 180° and 
360° refer to flou against the positive x-axis (i.e. coming from the open end). For the 
comparison with predictions, the experimental data were convened to this frame of 
reference. 


7.1. SPRE Test Case: 

Moderate Reynolds Number, Moderate Valensi Number 

7.1.1. Prediction of Laminar Oscillating Flow 

To demonstrate the capability of the numerical scheme adapted here, laminar oscillating 
flow, in a finite pipe was computed for Va * 80. The resulting velocity profiles in the axial 
center of the tube were compared with analytical results from the Uchida analysis. As can 
be seen from Fig 7.1, the agreement is excellent. The computations were made with the 
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^ 23 by 23 grid employed later to .be HRN k-e model. Alto 12 rime steps per cycle 
were used in connection with » Crank-Nicholson time imegnoon scheme. 

7.1.2. Predictions with the High-Reynolds Number k-e Model 

For preliminary studies. 3 cycles of oscillatory flow *> Re ma* - 1 1 700. Va -TO and 
L/D * 60 have been computed. Even though is known tha, the HRN tom of 
model (i) cannot predict transition and (it) is no. well suited to pipe flows a. moderate Re 
numbers, a computation seemed worfhwhile. First, the computation^ scheme could be 
checked widt the much simpler HRN formulation, and secondly, this computaoon proves 
, useful limiting estimate which can serve as an initial guess to more decried and exact 

computations. 

-The grid used here was a coarse grid of 23 by 23 grid points where the points in the 
nd ial direction were densely packed near the wall. Thus, to firs, internal grid potm was 
placed either in tire viscous sublayer or close to it throughout to cycle. The gnd was fine 
enough for .has test The firs, two cycles were computed widt 12 time steps per cycle, e 
dtird with 24. All computations were carried ou, wid, a fixed Crenk-Nicholson time 
integration scheme. The toee cycles come very Cose to. he periodic steady state as can be 
seen from to performance of to friction factor (Fig 7.2). A smaller time step does no, 
have a dramatic impact on the results of to primitive variables u and v. Thts suggests that 
to time stepping procedure is well chosen. A smaller time step, however, changes to 
values of k. e, more significantly than u and v. 

Since die HRN model does no, have the ability to predict transition, to turbulence 
mode, must be switched off "by hand" during the computations. In this study our emptncal 
cansttion criterion was simply to switch off the solution of to k and e equation whenever 
to Re number fell below 2300 and to switch i, on otherwise. For a Re number of less than 
2300. the laminar flow was anticipated, and the arrays of k. c and P, were set to aero. As 
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can be seen in the graphs below, the sudden disappearance of p. t at relative cycle times of 
0, 0.5, 1 etc. illustrate how unrealistic the transition process is modeled this way. 

It proved to be extremely difficult to keep the scheme stable at the point where Re « 0 
(reversal of the mean flow direction). Even after many trial runs it remained somewhat 
unpredictable whether the computation would "survive" the next mean flow reversal. This 
pointed out the need for optimizing the initial guesses from which a computation at a given 
time step was started and/or for another time stepping procedure. 

So far we used as initial guesses: 
p(i j) * 0 

u(ij) * (1-fgp) u 0 idCiJ) + fgp u m 

v(ij) = 0 

k(i j) = j* Tl*^um^ for Re > 2300 
e(i j) * 0.05 k^/v 

where fgp is a first guess parameter which varies continuously between zero and one (1 
for large time steps and few time steps per cycle, 0 for small time steps and many time 
steps per cycle). 

Also, reliable criteria for nonlinear convergence were determined. Define Omax the 
maximum absolute scaled value of the error in the mass conservation for any one control 
volume. Then: 


Omax < 


relax(u)*relax(v)*relax(p) 

ll*ml 


dOmax < 0 
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change in overall kinetic energy in domain ^ ^ 

relax(k)*relax(E)*relax(u t ) 

was required for at least 3 consecutive iterations. 

The boundary conditions for the turttulent kinetic energy and the turbulent dissipation 
„,e were 6% turbulence intensity and C - 0.05 in «,. (6.2). A velocity boundary c-dtocn 
of sinusoidally-varying slug flow into the domain was assumed. 

In the following, the results of the computations are presented. A comparison with the 
measured data of Simon, Seume and Friedman (1989) is given later together with the 
results of the low Reynolds number turbulence model. 

Results. The results of the computations are shown in Figs. 1.2 and 7.5 Figure 7.2 
shows the computed friction factor as a function of the crank angle and compares it with the 
faction factor cotrelations for steady situations at the respective Re number. 

Discussion. Compared to the laminar flow computation hardly any flow reversal is 
predicted This is good news for the boundary condition treatment of the outflow 
boundarv : Since at an outflow boundary outflow is assumed, complete upwinding is used 
there to eliminate the need to know the bounda^ conditions there. However, if there were 
to be inflow at an outflow boundary, information about the outflow boundary would be 
needed a priori. And this information is not normally available. 

AS can be seen from Fig.7.5. the gross charac.ens.ics of die flow cm be obtained fa. 
a HRN computation. But the near wall velocity and thus the friction factor predictions do 
not follow the experimental data well . This should no. be a surprise, because the validity 
of the universal law of the wall, on whtch the HRN relies, is very questionable for tire 

situation investigated. 
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Fig. 7.1 : Prediction of lami- 
nar oscillatory flow, Va * 80. 
Velocities correspond to a 
maximum pressure gradient 
which would yield a velocity 
of 1000 in steady flow 

: Uchida analysis. 

dot symbols : numerical 
computation. 


Fig. 7.2: Comparison of the 
friction coefficient prediction 
by the high-Reynolds number 
model with friction coefficient 
from turbulent and laminar 
steady state correlations. 



cronk angle 
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7.1.3. Prediction with the Low Reynolds Number k-E Model 

7. 1.3.1. Development of Computations 

For the SPRE test case (Remax = 1 1 700, Va = 80 and L/D = 60), a fluid flow 
computation was made using the Lam Bremhorst form of the low Re number k-£ 
turbulence model. For the first cycle computed with the LRN model, the computations at 
each time step were started from the results obtained with the HRN model described above. 
Except in the case where the number of time steps per cycle was 120, the computations of 
the later cycles were started from the corresponding results of the previous cycle. 

At first, four additional cycles were computed, cycles # 4, 5, 6 and 7, each with 24 
time steps per cycle. In cycle #7 the periodic steady state was reached and the results 
practically did not change any more. As inflow boundary conditions a TI of 5 % for the 
kinetic energy and a £ of 1/(TI Re® -875) [cf. gq. (6.20)] for the dissipation rate were used 

Later on, two more cycles (#8 and #9) were computed using measured velocity 
fluctuation data at the open end as the actual inflow boundary condition for k. This was 
done to eliminate the uncertainty associated with the assumed boundary condition in cycles 
#4 to #7. 

Finally, the grid independence of the results obtained was verified. 3/2 cycles were 
computed with 120 time steps per cycle and a grid of 33 by 51 grid points, and 3/2 cycles 
with 24 time steps and a larger grid of 64 by 91 grid points. The computation with 120 time 
steps per cycle was started at 0° crank angle with the results from cycle #9. The subsequent 
computations were started with the converged result of the previous time step. For the large 
grid computation, the results of cycle #9 were transposed to the finer grid, and the 
computations w ere started from the corresponding results of cycle #9. As will be shown 
later, the results of both tests were practically identical to cycle #9. The prediction of 
transition was not altered by the use of a finer gnd nor by the use of more time steps. This 
is an important finding since in Chapter 6 some grid dependence of the transition prediction 
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was found. The different behavior of the LRN model might be explained by the fact that in 
oscillating flow, due to history effects, the flow before transition is never totally 
undisturbed as assumed in the quasi-steady transition tests. From this we conclude that 
quantitatively trustful results may be, in principle, obtainable for oscillating flow from the 
turbulence model used. 

7. 1.3.2. A New Time Integration Scheme 

For the computation of the fourth cycle, the same time integration scheme as for the 
HRN model was used: a fixed Crank-Nicholson scheme. However, it was found that 
convergence frequently was not reached after as many as 200 nonlinear iterations. Then, if 
a solution of the next step was tried (even though the previous step was not fully 
converged), this solution nearly always diverged. Sufficient and dynamic under-relaxation 
and very many nonlinear iterations were the key to convergence at each time step. 

As initial guess, the values of the previous cycle were used at each time step of the fifth 
cycle. Initially it was assumed that only a few iterations would be necessary at each time 
step to reach convergence. But this hope did not come true. Again, many iterations needed 
to be done to prevent the subsequent time steps from diverging. 

A thorough review of the time integration scheme used so far revealed the underlying 
problem. In the LRN model computation 33*51 grid points w'ere used, most of them 
placed very close to the wall. As the numerical grid becomes finer and finer, the Crank- 
Nicholson formulation may become physically unstable. The fine grid near the wall 
practically assures that the the Crank-Nicholson scheme would become unstable unless 
many more time steps are taken. One possible solution would be to use a scheme between 
the fully implicit and the Crank-Nicholson scheme which is stable. This, however, would 
imply that the time integration scheme for the total domain would be based on the most 
unfavorable conditions for a Crank-Nicholson scheme in it, namely on the conditions near 
the wall. Near the wall, the effect of mass inertia is relatively small, whereas in the center 
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that effect is significant. Wherever inertia plays an important role, time history effects 
become important which are - in an always changing situation - better represented by the 
Crank-Nicholson scheme. Therefore, it is not as important to have a Crank-Nicholson 
scheme for the near wall region as it is for the center. This consideration bore out a time 
integration scheme which picks a time integration factor for each control volume and for 
each interface individually, based on the local conditions. For the situation under 
investigation here, the time integration factor varied from almost 1.0 near the wall to around 
0.95 in the center. Surely, for time integration factors that close to 1.0, a simple fully 
implicit method would have provided almost the identical results. This scheme proved to be 
stable at all times and was used for all cycles from #5 on. 

7. 1.3. 3. Results 

Comparison with measured uVms data. Assuming isotropic turbulence, the 
computed turbulent kinetic energy can be transformed to the rms axial velocity fluctuations. 

u' ms = V273T (7-2) 

Figure 7.3 shows a comparison of the measured velocity fluctuations with the 
computations at an axial location of x/D = 44 and 3 radial locations, centerline, intermediate 
and near the wall. The three computational curves plotted show the influence of the grid 
size in time and space. The measured T1 at the inflow is used in either case. However, even 
with a flat T1 of 59c (cycle #7) the computed curve looks alike and is not shown. The 
significant rise of u’ m s at circa 230° coincides in the experiment and the computations 
closer to the wall. However, the computations do not forecast a rise in the center. Also, the 
predicted decrease in turbulent fluctuations occurs over a much longer period than 
measured. It is believed that here the computation and the experiment show two complete!) 
different mechanisms of transition and relaminarization. The sudden decrease of u rms at 
300 c during the decelerating phase at a relatively high Re number indicates that the 
measured fluctuations between 230 c and 300 r correspond to a “turbulent slug” being sucked 
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in from the open end (cf. Seume, 1988). Without this slug, the flow at this location would 
remain laminar-like. On the other hand, the computations do not “see” that slug and rather 
describe an ordinary transition to turbulence at a too low, but computationally high enough 
Re number. This hypothesis is supported by Chapter 6.4., where it was found that in 
accelerated flow the turbulence model predicts transition at too low Re numbers. The 
computational transition predicts boundary layer instabilities and happens first near the 
walls. The high fluctuations closer to the wall will be transported to the center, but will 
reach the centerline only further downstream. Therefore, no rise of u’nns is shown at the 
centerline in the computations at this axial location. This explanation is supported by the 
computed fluctuations rising at about 90°. Here the axial location is further downstream of 
the inflow and the numerically predicted fluctuations spread over the entire cross section. In 
this case, experiment and computation “see” the same thing, i.e. ordinary transition to 
turbulence. Again, the turbulence model predicts transition at a too low Re number, given 
the rapid acceleration. The third rise in u’nns just after flow reversal at 180° is predicted 
very faithfully by the computation. It is believed that this rise is due to fluid that has 
become turbulent just after passing the probe location at x/D=44. Just after flow reversal, 
the fluctuations have not died down yet and revisit the probe location. 

In the following, we try to shed light on the question w-hy the turbulence model does 
not predict the experimentally observed turbulent slug, even though the correct inflow 
boundary condition is used for the turbulent kinetic energy. Fig. 7.4a shows the measured 
u’nns at inflow, the theoretical mean flow for this flow situation and the ratio of the two 
quantities u’ m s/um(t) which is equivalent to the TI at the inflow. The u’nns values were 
actually measured at the open end, but it is assumed that the inflow conditions are the same 
for both ends of the tube. The given curves can be repeated for crank angeles 180° to 360 c . 
Figure 7.4b shows axial profile of k at the centerline at different crank angles during the 
period of inflow from the open end. The values of k at x/D=60 correspond to the measured 
u’nns values of Fig.7.4a. It can be seen that the k values vary sharply between x/D=60 and 
x/D=58. At x/D=50, virtually all information about the inflow boundary condition is lost 
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Fig. 7.3: Comparison of computed velocity fluctuations with experiment a. 
open end). Top: r D-0.48; center: 0.433; bottom: centerline. 
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linear scale 



Seume and Friedman; computed with high-Re number model, dots indicate 


measured values;. LRN, 24 time steps per cycle, 33 by 51 grid, no experimental 

77 120 time steps per cycle, 33 by 51 grid, experimental Tl; 

LRX, 24 time steps per cycle. 64 by 91 grid, experimental Tl. Crank angles below ISO 0 
indicate a mean flow from the "piston end" of the experimental rig. crank angles betw een 
180 s and 360° indicate a mean flow from the open end. Vertical line indicates the centerline. 
To the right of the centerline , the velocity is plotted vs the radius linearly; to the left, it is 
planed against the logarithmic wall distance. The individual profiles are vertically shifted. 
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origins of shifted scale 


or dissipated. We suspect that inflow boundary condition for e provides unrealistically high 
values. High values of £ lead to a rapid dissipation of k just after the inflow and imply that 
in the computations no “turbulent slugs” could travel any significant distance downstream 
from the inlet. This suggests that further work is required to realistically specify e at the 
inflow and that £(Re) of eq. (6.20) should be modified. 

Comparison with measured velocity profiles. Figure 7.5 shows a comparison 
of the computed and measured velocity profiles at an axial location x/D = 44. It can be seen 
that the HRN and the LRN models give similar results in the center of the tube, but differ 
considerably towards the wall. In general, the HRN model predicts too high velocities near 
the wall. The results of three LRN computations are shown: 

(i) no experimental Tl, 24 time steps per cycle, 33 by 51 grid 

(ii) experimental TI, 120 time steps per cycle, 33 by 51 grid 

(iii) experimental Tl, 24 time steps per cycle, 64 by 91 grid 

The three LRN computations are practically indistinguishable. This verifies the grid 
independence and shows that the measured TI does not influence the predictions at this 
axial location. Tabic 7.1 tries to evaluate and to classify the results shown in Fig. 7. 5. 


Table 7.1: Evaluation of Figure 7.5 



Agreement with experiment 

Agreement with experiment 

Crank Angle HRN 

LRN 

Crank Angle 

HRN 

LRN 

30° 

satisf. 

good 

210° 

satisf. 

good 

60° 

fair 

satisf. 

240° 

fair 

satisf. 

90 

fair 

satisf. 

270° 

fair 

satisf. 

120° 

fair 

satisf. 

300° 

poor 

fair 

150° 

satisf. 

good 

330' 

fair 

fair 

180° 

good 

good 

360 

good 

good 


Scale: good - satisfactory * fair - poor 
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From this evaluation, a number of questions arise: 

1) Why is there good agreement for the LRN case at 150° and only fair at 330° ? 

2) Why is there satisfactory agreement for the LRN case at 120° and fair at 300° ? 

3) Why is there only little improvement at 60°, 90°, 300° and 330° from the HRN to the 
LRN model ? 

The following discussion shall assess these questions. From Fig. 7.5 and from Table 
7.1 it is evident that the results do not show symmetric deviations from the experimental 
results, i.e. the deviations at 150° differ significantly from the ones at 330° etc. For the data 
pair 607240° the predicted data is generally "too turbulent". However, this trend is more 
pronounced at 60°. The centerline velocities are underpredicted at 60° and on target at 240°. 
Near the wall, the absolute value of the velocity is overpredicted. An explanation for the 
deviation at 60° could be that the strong acceleration keeps the flow- longer laminar than the 
turbulence model can predict. This is in line with the test of the model for constant 
acceleration flows (see Chapter 6.4). From Fig. 7.3 it becomes clear that the better 
agreement at 240° comes from separated turbulent eddies being sucked in from the nozzle at 
the beginning of the second half of the cycle. The action of these eddies would increase the 
turbulent kinetic energy and counter the effect of the acceleration. Therefore, the measured 
data at 240° are "more turbulent" than those at 60°. This explanation also holds for the data 
pair 907270°. However, it appears that at 90° there is a very strong overprediction of the 
absolute values of the velocities between y/D=7xl0* 3 and 10 1 , whereas at 270° the 
prediction follows the experiment much better. Yet, at 270° all experimental data lie below 
or on the predicted curve. This indicates that the mass balance was not satisfied in the 
experiment. If the experimental data of 270° is shifted to give the same mass flow as in the 
computation, a similar trend as in the case of 240° can be seen. For 120° and 150° crank 
angle, similar deviation patterns are found. First, all experimental data are above or on the 
predicted curve, again indicating experimental differences in maintaining the mass flow 
rate. Near the wall, the prediction is right on target, whereas in the center, the predicted 
data is too low. The data at 300° and on 330° shows an overprediction of the absolute value 
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of the velocity near the wall, followed by an underprediction towards the center and the 
con-ect values at the centerline. The experimental data clearly looks laminar-like. Based on 
the measured boundary conditions for u’mis and on Fig- 7.3 that is due to the action of the 
nozzle before the flow enters the tube. There, fluid from the quiescent outside is sucked in 
and accelerated. The flow into the tube is relatively little disturbed and may well develop 
like a laminar flow initially. As seen in Chapter 6.2, at an axial distance of 16 diameters 
from the inlet, the flow will be laminar-like at moderate Re numbers unless the TI is very 
high. In the equivalent cases of 120° and 150°, the axial distance from the inflow is x/D = 
44. By then, a spatial transition is very likely to have happened which could explain the 
different behavior between 120° and 300° or 150° and 330°. 

It is remarkable that the near wall velocities are generally well predicted. It is the near 
wall velocities which determine the computed Cf value. Therefore the Cf predictions can be 
regarded with confidence. 

Given the remaining uncertainties of the experimental results, the general agreement of 
the LRN predictions with the experiment is considered to be good. 

Law of the wall. As can be seen in Fig. 7.6, the predictions support the hypothesis 
that the universal law of the wall is not a good representation for the velocity profile near 
the wall or even throughout the cross section. This does not come as a surprise since the 
universal law of the wall has already been shown for steady flows to be not applicable at 
low Re numbers. 

However, it can be seen that, except at flow reversal, there exists a laminar sublayer up 
to a y+ of about 7. Beyond this value, a logarithmic relationship between u + and y + may be 
formed, but the slopes are neither identical to the universal value nor constant at all. 

Tu and Ramaprian (1981) argue for pulsatile flow that the velocity does not scale with 
the w all shear stress at the same instance of time. Since the wall shear stress and the mean 
velocity have a phase difference could one scale the velocity with the shear stress of the 
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corresponding phase angle? To answer this question conclusively for oscillating flow, a 
phase relationship between u and to should be established, similarly like in laminar flow. 

Friction coefficient. The friction coefficient for fully developed flow derived from 
the turbulent steady state correlation does not agree well for the accelerated part of the cycle 
with the computed friction factor, which predicts lower values ( Fig.7.7). The agreement 
get rather close in the decelerated phase. The predictions of the HRN and LRN models are 
significantly different. Given the superiority of the velocity predictions obtained with the 
HRN model, it can be claimed that the HRN model does not give realistic values for Cf. 

Entrance length effects. Entrance length effects are important for about one third of 
the length of the tube during most of the cycle (Fig. 7.8). The fact that Cf,* initially 
decreases below the fully developed limit can be explained by the laminar-like flow 
development downstream of the inlet. This effect is relatively pronounced because of the 
low (experimentally determined) TI. It seems appropriate that a locally averaged friction 
factor for this case takes account for the entrance length effects. 

Other quantities. Fig. 7.9 shows the computed pressure distribution throughout the 
cycle. Figures 7.10, to 7.12 show the time variation of the turbulent kinetic energy, the 
dissipation rate and the turbulent viscosity at x/D = 44. Figure 7.13 is a vector plot of the 
the velocity. Figures 7.14, 7.15 and 7.16 show the variation of k, £ and (i-j at different 
crank angles. 
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Fig. 7.6: Comparison of near wall prediction of the low-Re number mode! with the 
universal law of the wall. x!D=44. Data point SPRE: Re max = 1 1700, Va-80, L'D-60 







Fig. 7.7: Comparison of 
friction coefficient prediction 
at the outflow cross section: 
Low-Reynolds number model 
vs. friction coefficient from 
turbulent and laminar steady 
state correlation. 



crank angle 


Fig. 7.8 : Ratio of local to 
fully developed friction 
coefficient at various crank 
angles. 
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Fig. 7.9: Predicted axial 
pressure distribution. Re ma x- 
11700, Va - 80, L/D = 60. 









Fig. 7.10: LRS predicted turbulent kinetic energy profiles at xlD - 44. Data point 
SPRE: Re max = 11700, Va = 80, L/D = 60. The half profiles do not come together at the 
centerline because the xJD location is not in the axial center of the tube 





Fig. 7.11: LRN predicted turbulent dissipation rate profiles at x/D = 44. Data point 
SPRE: Re max = 11700. Va = 80, L/D = 60. The half profiles do not come together at the 
centerline because the x'D location is not in the axial center of the tube. 
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Fig. 7.12: LRN predicted turbulent viscosity ' profiles at x/D = 44. Data point SPRE 
Re max = 11700, Va = 80, LID = 60. The half profiles do not come together at the centerliru 
because the x D location is not in the axial center of the tube. 







E: Re„ax= 1.1 7x1 0 4 . 









Fig. 7.15: Normalized turbulent dissipation rate eVrtr/Um.max at var ‘ ous crank angles. Data point SFRE: Remax=l .17x10* 
Va=80, UD=60. Mean flow direction: 120°. 150° into plane; 180° flow reversal; 210°, 240°, 270° out of plane. 








7.2. Other Test Cases Computed 


There were four more test cases computed, case e, d, m, p. The lettering of the cases 
corresponds to the names of Seume’s experiments. According to Seume’s findings, case e 
lies in the “fully turbulent” region, a region where the maximum Re number is high enough 
to cause instabilities to significantly perturb the flow, but where the frequency is too high to 
allow the fluctuations to die down as the flow reverses directions and accelerates. In this 
case, the ability of the turbulence model to predict transition is secondary. Case p is the 
corresponding case on the laminar side. According to the experiments by Seume (1988) 
and Ohmi et al. (1982) and Iguchi et al. (1983), the flow in case p always remains 


Table 72: Test cases investigated 


Case 

R^max 

Va 

Sir 

1VD 

Ar 

Mamax 

Ka,max 

SPRE 

1 . 1 7e4 

80 

0.0274 

60 

1.22 

0.015 

9.36e5 

d 

1.32e5 

81.2 

0.0025 

60 

13.6 

0.17 

1.07e7 

e 

1.87e5 

230.3 

0.0049 

60 

6.8 

0.23 

4.31e7 

m 

2.39e4 

230.3 

0.0386 

68.5 

0.8 

0.03 

5.50e6 

p 

8.43e3 

233.1 

0.0548 

60 

0.3 

0.01 

1.95e6 


laminar. Here, the ability of the turbulence model to accurately represent transition is the 
primary' factor for accurate predictions. Case d and m are in the transitional regime where 
the flow is laminar-like during parts of the cycle, and turbulent-like during the rest. In 
particular, with case d we can test the influence of increasing the Re number from the 
SPRE case while keeping the frequency constant; with case m we can test the influence of 
increasing the Valency number w hile maintaining the order of magnitude of Re number. 
The maximum non-dimensional acceleration occurring in each of the test cases falls 
approximately in the range of accelerations investigated in Chapter 6.4. 
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Comparison with experiment. No measured velocity profiles are available for 
either of these cases. However, since the profiles in the SPRE case can be predicted fairly 
well, we presume that also in the cases considered here the predicted profiles will be 
realistic provided that we can model transition correctly. For case p the laminar computation 
can give some basis for comparing velocity profiles. Available from experiment are the rms 
fluctuations at different locations for cases e, d and p. Even though case m is cited in 
Seume (1988), the data files could not be found any more. Figure 7.17 shows a qualitative 
comparison of the measured rms velocity fluctuations with the computed data at x/D=16 
(near the drive end). While the agreement is a very close between experiment and 
computation in cases e and d, it is decidedly worse for case p. 

In case e and d, the computed fluctuations in the center show a more structured, wave- 
like behavior than their experimental counterparts. In case e, the experimental near wall 
profile show s a clear phase lag as compared to the computations. As discussed in Chapter 
3.3.6., the relaxation time of turbulence seems to play a role in this case. The 
computational fluctuations are in phase with the mean flow field because the present 
turbulence model requires an immediate turbulent stress response to a large scale shear. The 
measured phase lag in the near w all fluctuations will most likely lead to a phase lag in the 
friction coefficient w hich is larger than predicted. In case d, where the frequency is only 
about one third of that in case e, no such phase lag can be detected in the measurements. In 
case p, the assumption of a specific value of TI at the inflow influences the result 
significantly. For a short pipe this outcome can be expected in light of the findings of 
Chapter 6. Based on the measured Tl in the SPRE case, it is believed that the 0.5 % TI for 
case p is closer to the experimental conditions. The computations show a clearer up-and- 
down trend than the experiment. Also, disregarding the pitfalls of a quantitative comparison 
for a moment, the level of the fluctuations predicted seems to be higher than seen 
experimentally. However, especially for the case of 0.5^r Tl, the level of fluctuations 
remains very low throughout, even near the wall. 
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Cross stream transport of turbulence.The effect of cross stream transport of 
turbulence is manifested by a phase shift between the near wall fluctuations and the 
centerline fluctuations and can be seen from Figure 7.18. The three cases shown (e, m and 
p) have the same frequency but different Re max - The computed shift is greatest for case p 
and small for e. In case e, practically during the entire cycle time, turbulence is generated 
near the wall and transported toward the center. In case p, fluctuations will be transported 
towards the center during the decelerating phase, and toward the wall in the accelerating 
phase. Case m lies in between those two extremes. 

Friction coefficient for fully developed flow. Figure 7.19 shows the 
computed friction coefficients and compared them with the laminar and turbulent steady 
state correlations. It becomes clear that for cases d and e the turbulent steady state friction 
coefficient is an excellent representation, whereas the results for case m are similar like for 
the SPRE case where the Cf values depart markedly from the steady state correlation in the 
accelerating phase. For case p, the steady state correlation is very bad. This, however, has 
long been known from analyses, experiments and computations of laminar oscillatory 
flow's. This finding supports the Shemer's (1985) hypothesis that a similarity parameter 
like the Va number which describes the influence of the unsteadiness of the flow on the 
various flow parameters like Cf should be built using some kind of effective viscosity 
instead of the molecular value. It is the effective viscosity that connects the motions of the 
boundary layer with the core in the tube. If the effective viscosity is high throughout the 
cycle (case d and e), then the influence of the unsteadiness on the flow parameters is small, 
even though the ordinarily used Va number suggests a strong influence. Note that for cases 
SPRE and m the average effective viscosity is higher in the decelerating phase where the 
steady state correlation agrees much better. Fig 7.19d shows, in addition to the laminar and 
turbulent steady state correlations, the results of a computation of case p in which the 
turbulence model was turned off completely. It can be seen that a variation in the 
specification ofTl influences the near-outflow Cf predictions only marginally. During the 
accelerating phase, the “turbulent" computations follow exactly the “laminar’' values. In the 
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decelerating phase, the “turbulent” computed friction coefficients are clearly higher than the 
corresponding “laminar” values. 

Entrance length effects. Effects of the hydrodynamical entrance length are 
negligible for cases e, d and m (Fig. 7. 20). Given the proximity of the Remax numbers of 
cases SPRE and m, it is striking that the entrance length effects arc quite different. 
However, the TI used in case m is a flat 5 % throughout the cycle which is much higher 
than what is used in the SPRE case. As known from Chapter 6, a higher TI causes the flow 
to develop sooner as a turbulent flow which explains this apparent discrepancy. Case p 
with 5% TI looks much like case SPRE. How ever, with 0.5% TI, the behavior at 150° 
crank angle deviates considerably from a steady state entrance length behavior. Generally, 
for case p, entrance length effects do play a role for the given pipe length. 

Other quantities. Figures 7.21 to 7.36 show 3-D plots for the axial velocity, 
turbulent kinetic energy, turbulent dissipation rate and turbulent viscosity for the various 
cases. It may be noted here that for case p the velocity distributions obtained with the 
turbulence model switched on look very much like the laminar profiles. 
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Fig. 7.17 a: Comparison of 
computed velocity fluctuations 
with experiment. Data point e: 
Re max = 1.87x1 0 5 . Va-230. 
xiD-lt (near drive end). 
UD=60. 



crank angle 


Fig. 7.1 7b: Comparison of 
computed velocity fluctuations 
with experiment Data point d 
Re ma x = ] -32x10 s . Va = 81. 
xD = 16 (near drive end), 
L'D=60. 



crcnk angle 
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Fig. 7.17c: Comparison of 
computed velocity fluctuations 
with experiment. Data point p: 
Re ma x=8-43x]0 3 , Va-231 , 
xlD-lt (near drive end). 
L/D=60. 



crank angle 
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Data Point c: Re=1.87e5 Va=230.3 



Data Point m: Re=2.39e4 Va=230.3 



Data Point p: Re=8.43e3 Vo=231.0 



crank angle 


Fig. 7.18: Radial transport of turbulence for three different data points with the same \ a 
number but different Re max number. 
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Fig. 7.19a: Comparison of 
computed fully developed 
friction coefficient with steady- 
state correlations. Data point 
e:Re max =l. 87x10 s . Va=230. 
UD=60. 
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Fig. 7.19b: Comparison of 
computed fully developed 
friction coefficient with steady 
state correlations. Data point 
d Re max = 1.32x1 0 5 . Va = 81. 
LD-60. 
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Fig. 7.19c: Comparison of 
computed fully developed 
friction coefficient: with 

steady state correlations. Data 
point in : Fe/nax — 2.39x10^, 
Va =230, UD=685. 



Fig. 7.19d: Comparison of 
computed fully developed 
friction coefficient with steady 
state correlations. Data point 
p Re ma x=S 48xJ0 3 . Va=23l, 
L'D=60 
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Fig. 7.20a: Ratio of local to 
fully developed friction 
coefficient at various crank 
angles. Data point e : 
Re ma x-l -87 x10 s , Va=230, 
UD=60. 


Fig. 7.20b: Ratio of local to 
fully developed friction 
coefficient at various crank 
angles. Data point d Remax - 
1 32x10 s . Va=81 . L/D=60. 





Fig. 7.20c: Ratio of local to 
fully developed friction 
coefficient at various crank 
angles. Data point tn : 
Re ma x=2.39xl0 4 , Va=230, 
UD=685. 



0.0 10.0 20.0 30.0 40.0 50.0 60.0 70. 


Fig. 7.20d: Ratio of local to 
fully developed friction 
coefficient at various crank 
angles. Data point p : 
Re ma x = 8 43xl0 3 , Va-231 , 
UD=60, 77=5.0%. 
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Fig. 7.20e: Ratio of local to 
fully developed friction 
coefficient at various crank 
angles. Data point p: 
Re ma x=S.43x]0 3 , Va=231, 
L/D=60, Tl=03%. 









Fig. 7.21: Normalized velocity u/u m rTUU at various crank angles. Data point e: Re max =l .87xl(P, Va=230. LID =60. Mean flow 
direction: 1 20°. 1 50° into plane; 180° flow reversal; 210°, 240°, 270° out of plane. 









Fig. 7.23: Normalized turbulent dissipation rate EV re r/Um.man at various crank angles. Data point e: Re max =l .87x10 s , Va=230. 
IJD=60. Mean flow direction: 120°, 150° into plane; 180° flow reversal; 210°, 240°, 270° out of plane. 






Fig. 7.25: Normalized velocity u/u m>ma x ot various crank angles. Data point d: Re max = 1 32x10 s , Va=8l, UD=60. Mean flow 
direction: 120°, 150° into plane; 180° flow reversal; 210°, 240°. 270° out of plane. 
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Fig. 7.29: Normalized velocity uAi m>max at various crank angles. Data point m: Re max =2J9xI0 4 , Va-230. UD-685. Mean 
flow direction: 1 20°. 150° into plane; 180° flow reversal; 210°, 240°, 270° out of plane. 







Fig. 731 : Normalized turbulent dissipation rate ev^f/u^ma* at various crank angles. Data point m: Re max =239xl0 4 . Va=230. 
UD=683. Mean flow direction: 120°, 150° into plane; 1 80° flow reversal; 21 0°, 240°, 270° out of plane. 









Fig. 7J3c : Laminar computation. Normalized velocity u/u m max at various crank angles. Data point p: Re max =8.43xl0 3 , 
Va=23l, L/D-60. Mean flow direction: 120°, 1 50° into plane ; 180° flow reversal; 210°, 240°, 270° out of plane. 



Fig. 7.34: Normalized turbulent kinetic energy k/(0.5 Um.max) at various crank angles. Data point p: Rtmax-^ ^xlO 3 . 
Va=23I . UD-60. Top row: T 1=5.0%; bottom row: TI-0.5%. Mean flow direction : 120°, 150° into plane; 180° flow 

reversal. 













7.3. Conclusions 


1. ) Starting from the results of the computations with the HRN turbulence model, the 

periodic steady state for case SPRE was successfully computed with the LRN 
turbulence model. In this model, no empirical transition criterion whatsoever was used. 

2. ) The predictions obtained are grid independent. 

3. ) Transition triggered by turbulent slugs (as seen in the SPRE case and described by 

Seume, 1988) is not predicted at the axial location investigated. In the SPRE case, 
ordinary transition is faithfully predicted but at somewhat too low Re numbers. The 
transition predictions for the high Re ma x cases e and d compare very favorably with the 
experiment. The always lamiar case p is computed very laminar-like. 

4. ) Based on the transition performance observed here, it is not necessary to take measures 

to broaden the predicted transitional Re number range (cf. Chapter 6.3.). 

4. ) For case SPRE, the computed velocity profiles at x/D = 44 agree rather well with the 

experimental data and show clear improvements over the HRN computation. 

5. ) The universal law of the wall does not hold for oscillating flow. However, a viscous 

sublayer following u + =y + does exist at least up to y + =7. 

6. ) The friction coefficient predictions show that for the two cases where Re max was 

greater than 1 0 5 , the steady state correlation is appropriate, at least up to Va=230. For 
the two cases where 10 4 < Remix < 10 5 , the steady state correlation can be used for the 
decelerating periods of the cycle. Here, the friction coefficients of the accelerating pans 
of the cycle have yet to be correlated. Below Re m ax = 10 4 the steady state correlation 
definitely does not hold throughout the cycle, at least for Va=200. 
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7. ) Entrance length effects are not important for the high Re m ax cases, but are significant 

for all other cases. 

8. ) To judge the impact of the unsteadiness on the flow parameters, an effective Va number 

should be used, Vacff 5 0 )R^/v c ff . 

9. ) In contrast to the findings of Rodi and Scheurcr (1986) for the flat plate boundary layer, 

here the LRN k-£ model does not seem to have particular problems with adverse 
pressure-gradients and decelerations. Rather, the problems come with strong 
acceleration. 

10. ) The shortcomings of the LRN computations are threefold: 

• In accelerated flow, the turbulence model predicts transition at too low Re numbers. 

• The inflow' boundary condition for £ does not reflect reality and does not allow 
turbulent slugs to exist long enough compared with the experiment. 

• The observed phase shift between the mean flow and the fluctuations at high Re and 
Va numbers is not predicted by the turbulence model. 
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PART III: heat transfer and irreversibility analysis 


8. HEAT TRANSFER IN STEADY PIPE FLOW 

In order to validate the pan of the computer program responsible for the heat transfer 
problem, laminar and turbulent steady pipe flow was computed for constant heat flux and 
constant wall temperature boundary conditions. The criteria for the validation were the 
following questions: 

1 . ) How well does the code predict Nusselt numbers for fully developed flow? 

2. ) Is the thermal entrance region realistically predicted? 


The Nusselt number correlations used to compare the predictions against were 1 : 



constant heat rate 

constant wall temperarure 

laminar flow 

Nu = 4.364 

Nu = 3.658 

turbulent flow 

Nu = 0.22 Re° 8 Pr° 5 

Nu = 0.021 Re 0 - 8 Pr0 5 


Figure 8.1 show s the computed local Nusselt numbers for a thermally and 
hydrodynamicall) developing pipe flow (L/D=150). The laminar computations were 
obtained w ith a 21 by 25 coarse grid, the turbulent computations with a 33 by 51 grid. The 
computed Nusselt numbers for the laminar case approach the theoretical values exactly. The 
computed values for the turbulent case are a little too high (136 vs. 121 for constant wall 
temperature, 138 vs. 126 for constant heat flux) but are within an error margin small 
enough to be acceptable. Also, the constant heat rate problem yields a higher Nu number 


'\V. M Kays and M E Crawford in Comeclhe Heat and Mass Transfer, MacGraw-Hitl Book 
Co., New Yo'k, 1980. 
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than the constant wall temperature problem. Yet, the fact that both numbers are over- 
predicted seems to indicate that the turbulence model employed slightly over-predicts the 
wall heat transfer. The entrance region is resolved realistically. For the laminar flow, the 
computed Nu number decreases monotonically to the asymptotic limit of fully developed 
flow. In case of turbulent flow, the predicted local Nu number decreases at first below the 
fully developed flow limit because of a short laminar-like development of the flow. As 
soon as the flow undergoes a spatial transition to turbulence (cf. Fig. 5.2 for the local 
frictions coefficient), the Nu number increases to the fully developed flow value. 

Figure 8.2 shows three-dimensional views of the temperature development in the pipe 
for two different dimensionless temperatures. In Fig. 8.2a the non-dimensional temperature 
is simply T/T 0 , whereas in Fig 8.2b it is (T w - T)/(T W - Tb u ik)- Those plots are supposed 
to be a reference against which one can compare the development of the temperature 
profiles in oscillating flow. 

For the 1 .5” pipe of the University of Minnesota oscillating flow experiment and air. 
the following relationship between the Ec number and the Remax number can be 
established: 


Ec = 5.2-10-13 

A0 

where A0 is defined as (T* - T in )/T m . Here, the representative temperature difference 
(AT)ref entering the Eckert number is taken as T w - T 0 , and T 0 = Tj n . Clearly, when A0 
approaches a very small value, the Eckert number becomes large and viscous effects play a 
significant role in the energy equation considered. In the given study, the temperature at the 
inflow was constant at a value of 300, and in the constant wall temperature case, T w was 
set to 360. This choice ensured that viscous heating was negligible. This is in line with 
Seume and Simon (1986) w ho stated that viscous heating does not play a role in Stirling 
engine heat exchangers. Despite that, the viscous dissipation function was always included 


in the calculations. 
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Fig. 8.1: Computed local Nu numbers and steady state correlations for fully developed 
flow. 



Fig. 8.2: Dimensionless temperatures T/T 0 and (T w - T)/(T W - T, n ) for Re-5u09(>. 
Pr-].0 and T w =constant. 
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9. PREDICTIONS OF HEAT TRANSFER IN 
TURBULENT OSCILLATING FLOW 


Given the assumptions outlined in Chapter 2, the properties are all considered as being 
constant in this study. Therefore the heat transfer problem is decoupled from the fluid flow 
problem and linear. Having obtained realistic solutions for the fluid flow problem, we can 
then expect realistic solutions for the heat transfer problem, too. The computations shown 
in this chapter will possess all weaknesses of the turbulence model discussed above. Those 
shortcomings will not be discussed here. 

For a practical application, the calculation of Nusselt numbers for the test cases 
considered is probably the most important task. The most precious question to be answered 
by this study is whether the steady state Nu number correlations are applicable in case of 
turbulent flow. Also of interest is the question whether entrance length effects play a more 
significant role in heat transfer than in the fluid flow part. We will restrict our attention to a 
Dirichlet boundary condition at the wall and inflow, and Neumann boundary condition at 
the centerline and outflow . The alternative case of a Neumann boundary 1 condition at the 
wall could easily be obtained, too. However, in case of turbulent flow and Pr=l, the 
difference in boundary conditions leads to an only insignificant difference of resulting Nu 
numbers. To maintain similarity with the University of Minnesota oscillating flow’ test rig. 
the temperature at the inflow cross section was assumed to be the same for inflow from the 
drive and from the open end. 

9.1. Nusselt Number Calculations near the Outflow 

Figure 9.1 show s the computed Nu numbers near the outflow cross section. As is the 
case for the friction coefficients, cases SPRE and m as well as cases e and d correspond to 
each other. Case p stands out alone. In all cases, the magnitude of the Nu number 
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throughout the cycle is much closer to the turbulent steady state correlation than for the 
laminar correlation. 

_y 

Cases SPRE and m, where Re max is moderate, display a definite phase shift of about 
the same magnitude (=20° crank angle). In the accelerating phase, the Nu number is less 
then the turbulent steady state correlation. This is expected since the accelerated flow is 
“more laminar-like”, and less cross stream eddy transport takes place. It may be noted here 
that the corresponding friction factors at 15° and 30° crank angle are larger than the 
turbulent stead)' state correlation, whereas the corresponding Nu number are smaller. For 
the friction factor this can be explained by the velocity gradients alone, which are very steep 
near the wall because of the acceleration. This is so even if there is no eddy transport. The 
Nu number not influenced by this velocity gradient, but determined by the eddy transport. 
On the other hand, the Nu number in the deceleration phase is enhanced by the increased 
eddy transport. 

The Nu numbers for cases e and d follow very nicely the turbulent steady state 
correlation. Apparently, the eddy activity here is high enough to lock the boundary layer to 
the core of the flow, similarly as in stead) flow. 

The Nu number pattern for case p shows the flow like laminar oscillating flow during 
the acceleration phase and deviates from that towards a more turbulent Nu number in the 
decelerating phase. The peak of the Nu numbers is offset by circa +80° compared to the 
steady correlation peak. 

9.2. Local Nusselt Numbers 

Figure 9.2 shows the local Nusselt numbers as computed. For cases e and d, where 
Remax is relatively high, the thermal entry length is short and the Nusselt number of the 
thermally and hydrodynamically fully developed flow gives a good representation for the 
entire tube. For cases SPRE and m, the thermal entiy length affects a significant portion of 
the tube. In the SPRE case, where the experimentally detennined, low TI boundary 
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condition is used, the initially flow develops laminar-like, and the Nu number is less than 
the fully developed value for about one third of the pipe length. In case m, a flat 5% TI at 
the inflow was assumed which causes the hydrodynamical entry length to be shorter than 
in the SPRE case. Here, about one fourth of the length is permanently a thermal entry 
length. The latter two cases show a similar behavior of the local Nu number at 30° crank 
angle: the local Nu number here is significantly higher than the fully developed Nu number 
practically throughout the tube. This can be explained as a history effect of the flow. What 
is the near outflow cross section during one half of the cycle is the near inflow cross- 
section during the other. Near the inflow, the respective inflow boundary conditions are 
strongly felt, whereas away from it the flow is more determined by local conditions. 

Shortly before flow reversal the level of turbulence is very low near the entrance. Then, 
after flow reversal, the level of turbulence at the former near-entrance cross section builds 
up only slow ly. In contrast, in the center of the tube or at the near-outflow cross section 
the level of turbulence is still relatively high before flow reversal . Therefore, just after flow 
reversal, the level of turbulence there is higher than at the now outflow cross-section. 
Consequently, the Nu number is lowest at the outflow cross section. In case m, due to the 
high frequency, this history effect is still slightly present at 60° crank angle. Case p (5% TI) 
is different from the previous two cases in that this history effect affects most of the cycle. 
Only the curves from 120° on are practically free of this effect. Hence, the fully developed 
Nu number value is not a good spatial mean. 

9.3. Temperature Solutions 

Underlying the Nusselt number results above are temperature solutions. In Figures 9.3 
to 9.5 the temperature solutions are shown in two different ways: The temperarure is 
normalized simply by division by a reference temperature T c . For the SPRE case, a 
normalization like in steady pipe flow is used, G = (T* - T)/(T W - Ttj U u<j. The shown plots 
look similar like plots for steady state, but the normalization brakes down for flow reversal. 
Only cases SPRE. e and p are shown. 
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9.4. Conclusions 


1. ) For Re max numbers above 10 5 , the turbulent steady state Nu number correlation 

approximates well the computed instantaneous Nu number for the fully developed 
flow, at least up to Va=230. 

2. ) For 10 4 < Remax < 10 5 , the fully developed instantaneous Nu number can be related to 

the steady state correlation by a simple phase lag relation (to be developed). 

3. ) Below Re m ax * 10 4 , the instantaneous fully developed Nu number differs in phase and 

magnitude from the turbulent steady state correlation. 

4. ) Thermal entry length effects are negligible for Re m ax numbers above 10 5 . 

5. ) Below a Remax of 10 5 , the thermal entry length becomes appreciable and history effects 

begin to play a role. 
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Fig. 9.1a: Comparison of 
computed fully developed hlu 
number with steady state 
correlations. Data point 
SPRE: Re ma x=l .17xl0 4 , 
Va=80,UD=60. 



crank angle 


Fig. 9.1b: Comparison of 
computed fully developed Nu 
number with steady state 
correlations. Data point m: 
Re mai = 2.39xl0 4 . Va=230 . 
UD=685. 
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Fig. 9.1c: Comparison of 
computed fully developed Nu 
number with steady state 
correlations. Data point e: 
Re max = 1.87x10*. Va-230, 
UD=60. 
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Fig. 9. Id: Comparison of 
computed fully developed Nu 
number with steady state 
correlations. Data point d: 
Remax-1 .32x10*. Va=81 . 
LD=60 
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Fig. 9.1e: Comparison of 
computed fully developed Nu 
number with steady state 
correlations. Data point p: 
Re maj =8.43xl03, Va-231, 
UD=60, Tl-5.0%. 
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Fig. 9.1 f: Comparison of 
computed fully developed Nu 
number with steady state 
correlations. Data point p: 
Re max -8.43xl0^ , Va-231 , 
UD=60, Tl^OJVc. 
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Fig. 9.2a: Ratio of local to 
fully developed Nu number at 
various crank angles. Data 
point SPRE: Re max - 
1.17x10* , Va-80, UD=60. 


Fig. 9.2b: Ratio of local to 
fully developed Nu number at 
various crank angles. Data 
point m: Re max =2 .39x1 0* . 
Va=230, UD=685. 










Fig. 9.2c: Ratio of local to 
fully developed Nu number at 
various crank angles. Data 
point e: Re ma x = l 87x105, 
Va=230, UD-60. 


Fig. 9.2d: Ratio of local to 
fully developed Nu number at 
various crank angles. Data 
point d: Re ma x = 1 .32x1 0$ , 
Va=81. UD=60. 







Fig. 9.2e: Ratio of local to 
fully developed Nu number at 
various crank angles. Data 
point p: Re max =8-43xl0 3 , 
Va-23J, L/D-60, Tl=-5.0%. 


Fig. 9.2 f: Ratio of local to 
fully developed Nu number at 
various crank angles. Data 
point p: Re ma x = 8 -43xl0 3 , 
Va=23), LiD-60, Tl=OJ%. 
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Fig. 9.5: Normalized temperature TIT 0 at various crank angles. Data point p: Re m<u =8. 43x/0*, Va=231, UD=60, Pr=l .0. 
F.c—tQ, TJT 0 -].2, Tl=Q5. Mean flow direction: 120°, 1 50° into plane; 180° flow reversal; 210°. 240°, 270° out of plane. 



10. ENTROPY GENERATION IN STEADY PIPE FLOW 


10.1. Derivation of the Entropy Generation Term 

As outlined in Chapter 2, the benefit of the differential equation for the entropy lies in 
its entropy generation term. In heat exchanger design, techniques to reduce friction and to 
enhance heat transfer are usually in conflict with each other. Applying the second law of 
thermodynamics puts both irreversible processes on the same physical scale and allows to 
properly evaluate their impact on the overall performance. 

The differentia] equation for entropy can be written as 1 (Bejan, 1982) 


p § = - div( Y ) + s gen 


( 10 . 1 ) 


the energy equation is 

p Ml = -div( q ) + p<J> 
and the continuity equation is 


f^ + pdiv(u) = 0 


The first term on the right hand side of eq. (10.1) can be rewritten as 


( 10 . 2 ) 


(10.3) 


1 In ihe following, ihe quamiu 


Di 


denotes the substantial derivative. 
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(10.4) 


div( y ) = ~ y div( q ) + • grad(T) 


Replacing div( q ) with the energy equation, the equation for the entropy becomes 


pB^A.gradCn+l 


iLph_2p 


P ^r-TT 
K Dt Dt 


+ $ 


gen 


(10.5) 


Using Gibbs’ equation T ds = dh - (1/p )dp and Fourier’s law, eq. (10.5) transforms to 


\ (P at) (T)] 2 + Y ^ 


(10.6) 


For turbulent flow , the ensemble averaged equations are 


P § = ” d ’ V( T ) + ^* e " ’ P d ’ V( “ S ^ 


(10.7) 


— = -div( q ) + ^ + pC> + pe - p div( iT h ) + div( u p ) 


( 10 . 8 ) 


It is 


div( u s ) = s div( u ) + u • grad( s' ) (10.9) 

and similarly for h’ and p’. For constant density flows, the first term on the right hand 
side drops out. 

Combining equations (10.7) and (10.8), using (10.9) and applying Gibbs’ equation for 
the fluctuation terms, the entropy generation term for turbulent flow takes the final form 
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( 10 . 10 ) 


S,„ = ilg^(T>] J + = 4 > + =' 

T 1 ] 

The individual terms in eq. (10.10) reflect the irreversibilities due to heat conduction, 
mean flow across a finite pressure gradient and fluctuations across finite pressure 
gradients. All terms of eq. (10.10) are positive which is in accordance with the second law. 
It is noteworthy that the heat transfer term in eq. (10.10) is governed by the molecular , not 
the effective conductivity. 

Using the non-dimensional variables as introduced in Chapter 2 and defining 


T 



( 10 . 11 ) 


£ = 


p R 




( 10 . 12 ) 


T D 

c = 1° S 

Jgen ~ 2 8 en 


P re f ^m.max 


(10.13) 


the non-dimensional equation for the enttopy generation rale becomes: 


** en ~ PrEc (T* 


T ) 

in 


-L[grad(T)f + A ?;£ 


(10.14) 


where the factor 4 on the left hand side is due to the use of D as length scale for the rate 
of generated entropy. Here, the case Ec = 0 represents a singularity, and it is not clear a 
priory how small Ec must be in order to justify an omission of the two last terms. 
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To estimate how much entropy is generated in a numerical control volume, the 
generated entropy of eq. (10.13) will be multiplied with the dimensionless volume of the 
control volume. 


5 


gcn,num 



AV cv 
0.5 R 2 L 


(10.15) 


Note: The entropy production depends on the absolute temperature at which the 
irreversible processes occur. Therefore, all results shown are only valid for the absolute 
temperatures picked, or, more exactly, for the ratio of the absolute temperatures. 


10.2. Results for Steady Turbulent Flow 

Figure 10.1 a shows the total non-dimensional entropy generation rate for a steady 
turbulent pipe flow at Re = 50 000, assuming Tm = 300 and T w = 360 (like in the heat 
transfer problem). Most of the entropy production occurs near the wall and near the 
entrance cross section where the gradients are especially steep. 

Figure 10.1b shows the ratio of thermal entropy production (first term on RHS of eq. 
(10.10) ) to total entropy generation. Figure 10.1c shows the ratio of frictional entropy 
production to total entropy production, and Figure 10. Id shows the ratio of turbulent 
entTopy production to total entropy generation. Near the wall, where the most entropy is 
produced, the thermal production is dominant. Frictional production is negligible. 
Turbulent entropy production is significant towards the center of the tube and near the 
entrance. However, of all three contributors, thermal entropy production is largest. 
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Fig. 10.1: Entropy production rate for Re = 50000, Pr=1.0, T w lT 0 =l .2, T ir = T c . 

Top left: Normalized total entropy generation rate; top right: thermal fraction oj total 
generated entropy: bottom left: frictional fraction; bottom right: turbulent fraction. 
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11. ENTROPY GENERATION IN TURBULENT OSCILLATING 

FLOW 


For the different test cases, it is of interest for the designer of a heat exchanger to have 
the following two questions answered: 

1. ) Where does the significant portion of the irreversibility take place? 

2. ) Which process contributes most to the irreversibilities, heat transfer, friction of the mean 

flow or turbulent dissipation? 

As outlined in Chapter 10, the non-dimensional entropy depends on the absolute 
temperature chosen to non-dimensionalize it. Since there is no compelling reason to choose 
a particular temperature, its selection is somewhat arbitrary. This is a well known problem 
in exergy analyses, too. However, from exergy considerations it can be argued that the 
ambient temperature (if any clear defined T» exists) is the preferred choice. Hence, the 
results shown here are only exemplary for one typical case where Tc^T* and T w /Toc=1.2. 

Results and Discussion. Figures 11.1 to 11.3 show the total normalized entropy 
generation in the domain at different crank angels for cases e, SPRE and p. In each case, 
nearly all of the entropy is generated very close to the wall. The peak generation is very 
close to the entrance cross section. 

Figures 1 1 .4 to 1 1 .6 show the portion of entropy generation due to heat conduction for 
cases e, SPRE and p. It becomes clear that, overall, conduction is the main contributor to 
irreversibilities. This can be explained by the fact that the thermal irreversibilities depend on 
the gradient of T, whose radial component is zero at the centerline. The turbulent 
dissipation, in contrast, influences irreversibilities directly and is nonzero at the centerline. 

Irreversibility contributions of turbulent dissipation do become more significant towards 
the center and in the inflow region (Figs 1 1 .7 to 1 1.9). Especially in case e, turbulent 
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dissipation seems to be the main contributor^ %) for much of the domain. However, since 
also in this case by far the most overall entropy production takes place in a very thin layer 
near the wall, thermal entropy production is the largest contributor to irreversibilities. 

Conclusions. No generally valid conclusions can be drawn. For the cases considered 
here, thermal entropy production was largest. However, this constellation can change if the 
temperatures involved change, or if the Remax number is increased significantly. The above 
posed questions must be answered individually from case to case. This points at the need 
for reliable computer programs w ith w hom each case can be simulated separately. 
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PART IV: CLOSURE 


12. OVERALL ASSESSMENT 
12.1. Summary and Major Conclusions 

1. ) A literature review shows that presently no detailed numerical study into fully 

developed and developing turbulent oscillating flow and heat transfer is available. 

2. ) A control-volume based numerical algorithm suitable for solving the governing 

equations exactly and efficiently is developed. 

3. ) The k-e model in the Lam-Bremhorst form is identified as a suitable model for 

oscillating flow predictions. It is shown that the model has the capability to predict 
transition to turbulence in quasi-steady and accelerated pipe flow at least qualitatively 
correct. 

4. ) The oscillating flow predictions generally compare well with the experiment. This 

validates the choice of the k-e model for this study. 

5. ) With regard of the oscillating flow predictions, the major flaws of the k-e model are: 

• The e inflow boundary condition is questionable. 

• Transition in accelerated pipe flow is predicted too early. 

• The present eddy-viscosity concept implies an infinite stress-response time to shear. 

6. ) A modification to the e-equations is proposed in order to capture better the 

accelerarion/deceleration effects on transition and relaminarization. 

7. ) A complex valued turbulence viscosity is proposed. 



8. ) The universal law of the wall does not hold for oscillating flow. A viscous sublayer 

following u + =y + does exist at least up to a y+ of 7. 

9. ) For Remax numbers above 10 5 , the steady state correlations for the fully developed 

friction coefficient and Nu number apply well. For lower Rema* numbers, the 
departures become larger with decreasing Remax and increasing Va. 

10. ) To estimate the influence of unsteadiness, an effective Valensi number should be used, 
built upon the effective, not the laminar, viscosity. 

1 1 . ) An irreversibility analysis demonstrates that for the conditions chosen, heat conduction 
is the biggest contributor towards entropy production. 

12.2. Contributions of This Research 

To the best knowledge of the author of this work, the new and unique features 
contained herein are: 

1. ) The development of the locally adaptive time integration scheme for a nonlinear 

convection-diffusion situation. 

2. ) The application and documentation of the predictions of the Lam-Bremhorst form k-£ 

model for quasi-steady and accelerated fully developed and developing pipe flow. 

3. ) The oscillatory flow and heat transfer predictions. 

4. ) The proposed modification to the E-equation. 
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12.3. Suggestions for Further Research 

Four major poinis emerged for further investigation: 

1 . ) The inflow boundary condition for £ should be theoretically or, if possible, 

experimentally investigated in order to enable the prediction of “traveling turbulent 
slugs” downstream of the inflow. 

2. ) The concept of a complex valued turbulent viscosity should be pursued to capture the 

phase lag of the small scale motion with respect to the large scale motion in strongly 
unsteady flows. 

3. ) The proposed modification to the e-equation should be tested and scaled in order to 

yield better predictions for accelerated/decelerated flows. 

4. ) The data generated herein should be reduced to yield correlations for the friction factor 

and Nu number which are needed in Stirling engine performance codes. 
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APPENDIX 


A. Vector Quantities in Axisymmetric Coordinates 


Note: In this appendix, the following conventions apply: 
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Kinematic vector quantities in axisvmmetric coordinates! 
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div[ngrad(u)f] = 



Viscous dissipation function <t> and generation term G for turbulent kinetic energy: 


. System of Differential Equations Solved 
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C. Grid Generation 


The grid generation can be divided into the two parts of dividing the radial and axial 
directions into small stretches from which the control volumes will be built. The width of a 
control volume will be denoted as YV(J) for the radial direction and XU(I) for the axial 
direction. There are Ml grid lines dividing the radial direction and LI for the axial 
direction. J is the numbering of the grid points for the radial direction, I for the axial 
direction. The axial grid lines were equidistant in all computations. 

The objectives for the radial grid were: 

• very dense near the wall 

• neighboring control volumes should not vary too much in width 

• sufficient number of grid points also in the center region should be maintained 

Let us define x radial grid coordinate 

J - 2 

y e 

4 Ml - 2 

and Xi be an intermediate location where we switch from one grid form to another. 
There are many approaches thinkable. However, in axisymmetric coordinates, the 
following customary approach does not work well: 

YV(J) I a * +b * <J£ > 

YL l cx d + e % >X t 

The reason for this is that in axisymmetric coordinates, not enough grid lines can be 
placed in the region near the wall--which should be resolved finely-since the quantity x get 
close to 1 .0. even with a high exponent d. On the contrary, near the centerline, x 
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approaches zero, and ihe resulting variation in the width of neighboring control volumes is 
tremendous. 

A remedy is to specify the following grid: 


F(X) = 


YV(J) 

YL 


ax X < Xj 
a e bx + c x>Xj 


Which gives a linear grid between the centerline and Xi and an exponential grid 
between xi and the wall. The parameters b and xi must be specified, the others are fixed. 

At the wall, x = 1 =>. F(x) = 1. Therefore: 


c = 1 - a e b 

At xi. require that the grid function is continuous and differentiable. Thus 

— bTT 

e b + e (bXj-1] 

and a = a b e^ x i^ 

For b — » - °° =0 very dense near the wall; 

for b — * + => very dense near Xj- 

This form was used for the grid generation of the LRN computations. For the 
computations with the 33 by 51 grid, xi = 0.4 and b = - 5 were customarily used. For 
the computations with the 35 by 64 grid, xi = 0.2 and b = - 6 were used. Figure C.l 
shows a typically used grid. 


229 



grid of r910125 

y-dimension stretching by a factor of 120 



Figure C.I: A typical grid, 64 by 91 grid points. Mote that the radial dimension is 
stretched by a factor of 120. 
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